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1.  Introduction 

Scores  of  research  activities  have  been  conducted  in  recent  years  as  part  of  an  attempt  to 
define  the  propagation  and  attenuation  of  communications  signals  through  the 
atmosphere.  These  have,  to  a  certain  extent,  allowed  for  the  consideration  of  numerous 
factors,  including  multipath,  ducting,  diffraction,  terrain,  and  weather  phenomena. 
However,  while  it  is  the  case  that  RF  propagation  models  can  provide  a  characterization 
of  the  physical  aspects  of  these  signals,  future  projects  will  require  realistic  and  real-time 
representations  of  communications  signals  with  these  physical  aspects  incorporated.  The 
Air  Force  Flight  Test  Center  and  the  Navy  Air  Combat  Environment,  Test  and 
Engineering  Facility  both  require  the  capability  to  accurately  represent  the  “real  world” 
atmospheric  propagation  of  communications  signals  for  real-time  Installed  Systems  Test 
Facility  (ISFT)  applications. 

This  research  effort  is  directed  at  the  fulfillment  of  this  demand  through  the 
accomplishment  of  the  following  tasks:  Search  of  the  available  literature  and  models  of 
these  phenomena;  identify  those  numerically  efficient  algorithms  and  models  with 
application  to  the  installed  systems  test  environment;  and  lastly,  specify  and  prototype  a 
software  module  for  incorporation  into  the  existing  communication  simulator.  The  effort 
will  maximally  utilize  existing  research  and  modeling  and  further  extend  those  efforts  to 
the  ISTF  application. 

Specifically  this  task  will  develop  methods  for  rapid  efficient  methods,  algorithms  and 
software  to  signal  attenuation  for  the  Joint  Communications  Simulator  (JCS).  The  JCS  is 
a  real  time  simulator  and  thus  the  delivered  software  must  provide  propagation 
predictions  in  real  time  for  all  relevant  simulation  players.  The  Test  and  Evaluation 
scenarios  consist  of  a  400  x  400  square  km  battle  space  area  in  which  the  system  under 
test  is  immersed.  Every  simulation  player  is  moved  on  1/20  second  time  steps  producing 
a  typical  aircraft  movement  of  10  meters  per  time  step.  Since  there  are  up  to  400  mobile 
platforms  and  2000  emitters  in  the  scenario,  this  gives  a  potential  required  execution  time 
for  a  single  propagation  prediction  of  no  less  than  1/80,000  seconds.  The  calculation  time 
constraint  includes  not  only  determining  the  physical  interactions,  but  also  extraction  of 
terrain  data  and  other  environmental  data.  Obviously,  if  possible,  only  those  propagation 
predictions  in  which  the  loss  displays  a  large  change  over  a  typical  time  step  should  be 
evaluated.  Thus,  the  terrain  resolution  and  other  first  order  environmental  factors  will 
drive  determining  the  re-computation  time  for  the  simulation  assets. 

The  priorities  for  this  research  were  coordinated  with  the  JCS  program  manager  and  the 
developers.  Their  foremost  concern  is  giving  JCS  the  capability  to  model  terrain  effects. 
Initially  the  priority  is  for  first  order  terrain  effects  caused  by  terrain  blockage.  Once  the 
first  order  attenuation  effects  have  been  integrated,  terrain  multi-path  is  to  be  considered. 

Because  the  JCS  currently  only  computes  free  space  loss  without  considering  terrain 
blockage  or  reflection,  the  first  phase  of  this  research  effort  shall  explore  computational 
methods  for  Line  of  Sight  (LOS)  interference  and  Beyond  Line  of  Sight  (BLOS) 
diffraction  effects.  Further  the  JCS  has  only  96  signal  generators  available  to  model 
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multi-path  signals,  thus  only  the  direct  and  a  representative  composite  reflected  ray  will 
be  used  for  the  simulated  communications  link.  In  addition  to  improving  LOS  and  BLOS 
methods,  optimal  methods  to  extract  terrain  data  while  determining  propagation  relevant 
geometrical  properties  shall  be  investigated.  During  this  phase  we  intend  to  improve  upon 
the  accuracy  and  run  time  speed  of  the  methods  commonly  found  in  the  Joint  Spectrum 
Center’s  (JSC)  Terrain  Integrated  Rough  Earth  Model  (TIREM)  and  the  Institute  of 
Telecommunications  Service’s  (ITS)  Irregular  Terrain  Model  (ITM),  also  known  as 
Longley-Rice.  Parabolic  Equation  methods,  unless  implemented  on  signal  processing 
hardware,  have  been  deemed  to  be  much  too  slow  to  be  used  for  the  JCS. 

The  second  phase  of  the  project  seeks  to  extend  the  methods  from  phase  1  to  be  able  to 
model  multi-path.  This  will  expand  upon  the  single  reflected  ray  method  to  give  a  more 
statistical  representation  of  the  multi-path.  Although  many  methods  are  currently  applied 
to  model  multi-path  in  signal  simulations,  the  computation  constraints  will  force  the 
multi-path  model  to  be  very  fast  and  use  limited  computer  resources.  Further  the  multi- 
path  methods  will  not  be  able  to  utilize  the  JSC  signal  processing  hardware. 

In  short,  Remcom  shall  develop  efficient  methods  to  model  terrain  reflection  and 
diffraction  including  multi-path.  Operational  software  based  on  the  model  will  be 
developed  in  order  to  provide  JCS  the  capability  of  running  faster  than  real  time 
propagation  predictions  required  by  the  test  and  evaluation  scenarios. 

2.  Status  of  effort 

This  model  will  not  only  be  used  by  the  ISTF  JCS,  but  also  by  the  Air  Force  Intelligence 
Warfare  Activity  (AFIWC).  Because  the  AFIWC  requires  a  similar  model  with  an 
extended  frequency  range  to  500  KHz,  AFIWC  is  funding  the  development  of  the  LF  and 
HF  ground  wave  methods. 

This  report  summarizes  the  accomplishments  toward  the  goals  of  the  first  phase.  Since 
the  last  status  report  in  September,  2003,  we  have  continued  testing  our  methods  and 
evaluating  their  performance  in  preparation  for  integration  with  the  terrain  extraction,  and 
transition  to  the  JCS  architecture.  The  project  plan  is  found  in  Section  6  and  the  spending 
curves  found  in  Section  7. 

The  prototype  of  Remcom’s  new  propagation  model  incorporating  first  order  terrain 
effects  has  been  completed.  This  model,  referred  to  as  the  Terrain  Multipath 
Electromagnetic  Model  (TMEM),  is  described  in  detail  in  the  next  section. 

3.  Terrain  Multipath  Electromagnetic  Model 

Existing  propagation  models  force  a  tradeoff  between  run-time  and  accuracy.  Models 
based  on  the  parabolic  equation  (PE)  such  as  APM  and  VTRPE  provide  good  results,  yet 
are  computationally  intensive.  Other  methods  are  quicker,  such  as  those  used  by  TIREM 
and  ITM,  but  results  less  realistic. 
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Remcom  has  developed  the  Terrain  Multipath  Electromagnetic  Model  (TMEM)  to  fill  the 
need  for  a  propagation  model  that  is  fast  enough  for  simulation,  yet  provides  realistic 
results. 

3.1.  Definitions 

Path  Loss  is  defined  as  the  log  of  the  ratio  of  the  power  actually  received  at  a  specified 
antenna  to  that  transmitted. 


L  =  -101og10^— j 

The  received  power  Pr  may  be  expressed  (MKSA  units)  as  the  product  of  the  power 
density  of  the  electric  field  at  the  receiver  ( watts/m 2)  and  the  theoretical  maximum 

aperture  area  available  to  the  emitter  which  for  a  Herzian  dipole  is 
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and  the  constant  radiated  power  at  the  transmitter  is  the  product  of  the  corresponding 
power  density  and  the  full  surface  area  of  the  sphere  at  radius  r,  the  line-of-sight  distance 
between  transmitter  and  receiver. 
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If  received  and  transmitted  fields  are  equal,  the  ratio  of  fields  in  the  second  bracketed 
expression  is  unity  and  the  path  loss  reduces  to  the  first  bracketed  term  which  defines  the 
free  space  loss.  The  ratio  of  received  to  transmitted  fields  defines  the  attenuation  factor, 
and  it  is  the  principle  objective  of  terrestrial  propagation  to  model  the  log  of  the 
attenuation  factor  in  dB  to  obtain  the  loss  over  free  space  or  excess  loss. 

In  this  proposal  the  diffractive  effects  defining  the  excess  loss  are  modeled  as  products  of 
two  elementary  forms  including  that  for  diffraction  over  irregular  terrain  and  associated 
with  diffraction  by  a  perfectly  absorbing  knife  edge  and  that  for  smooth  earth  propagation 
and  modeled  according  to  the  numerically  accessible  methods  of  Norton  [10].  The 
profiles  will  be  subdivided  according  to  the  geodesic  geometry  of  the  terrain  and  each 
subdivision  will  be  categorized  according  to  which  of  the  two  elementary  diffraction 
models  it  best  conforms.  The  total  loss  will  be  defined  as  the  product  of  partial 
attenuation  factors  resulting  from  these  elementary  assignments  and  of  a  functional  form 
intended  to  achieve  continuity  as  the  various  partial  diffraction  sectors  are  traversed.  The 
following  section  will  display  the  essential  features  of  the  terrain  subdivisions  according 
to  elementary  diffractive  forms  and  will  include  examples  with  smooth  earth  terrain  and 
an  example  consisting  exclusively  of  irregular  terrain  features. 

3.2.  Comprehensive  Perspective  (Overview)  of  Partial 
Diffraction  Factors  Relevant  to  Terrestrial  Propagation 

The  elementary  diffraction  forms  discussed  in  the  previous  section  may  be  identified  by 
the  presence  of  characteristic  features  of  the  geodesic  geometry  inherent  in  the  terrain 
profile.  As  applied  to  geometrical  optics,  geodesic  paths  are  the  ray  paths  of  shortest 
distance  from  a  fixed  transmitter  to  a  receiver  of  arbitrary  range  and  height  along  the 
terrain  profile.  The  diffraction  centers  used  to  compute  partial  attenuation  factors  for  both 
irregular  and  smooth  terrain  may  be  identified  as  the  vertices  of  the  discrete  portion  of  the 
geodesic  path  or  as  end  points  of  the  tangent  arcs  accompanying  terrestrial  geodesic  paths 
over  uniformly  elevated  terrain.  The  (discrete)  diffraction  centers  are  accompanied  by 
optical  accessibility  sectors  within  which  unobstructed  ray  paths  exist  between  the 
diffraction  center  and  a  field  point  confined  to  the  sector.  The  angular  bounds  of  the 
sectors  may  be  defined  such  that  they  subdivide  the  vertical  propagation  plane  into 
regions  dominated  by  the  defining  diffraction  center  which  is  located  at  the  intersection 
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of  the  ray  paths  defining  its  upper  and  lower  limits.  An  example  of  the  subdivided  radial 
propagation  plane  for  a  path  profile  consisting  exclusively  of  irregular  terrain  appears  in 
Figure  1. 


RANGE  COORDINATE  (KM) 

Figure  1.  Subdivide  Radial  Propagation  Plane  for  Irregular  Terrain 

3.2.1.  Discussion  of  irregular  terrain  diffraction  sectors 

The  terrain  profile  appearing  in  Figure  1,  consisting  of  a  random  assortment  of  knife 
edges,  will  be  used  to  illustrate  the  diffraction  sectors  described  above  by  examining  the 
functional  behavior  of  the  total  path  loss  in  a  vertical  excursion  along  the  line  segments 
appearing  near  the  right  margin  of  the  figure.  The  total  attenuation  due  to  the  knife  edges 
is  modeled  as  the  product  of  single  knife  edge  attenuations  in  which  the  effective  source 
(transmitter)  geometry  is  provided  by  the  most  immediate  geodesic  relation  opposite  the 
direction  of  propagation.  The  total  loss  is  then  the  sum  of  the  (partial)  losses 
corresponding  to  each  term  in  the  product  of  attenuations. 

A  set  of  vertically  arrayed  symbols  has  been  introduced  in  the  figure  for  each  diffraction 
sector  to  indicate  the  evolution  of  its  partial  excess  loss  contribution  as  a  function  of 
height.  The  diffraction  sectors  in  Figure  1  have  been  labeled  with  the  same  Arabic 
numeral  as  that  used  to  label  the  associated  knife  edge  (diffraction  center).  The  geodesic 
path  to  an  arbitrary  field  point  is  a  simple  extension  of  the  path  to  the  diffraction  center 
along  a  line  segment  from  the  diffraction  vertex  to  the  field  point.  The  specific  partial 
attenuation  contributed  by  the  sector  depends  upon  the  location  of  the  field  point  relative 
to  that  sector.  Thus  the  boundary  between  sectors  0  and  1  is  marked  by  a  solid  line 
parallel  to  the  shadow  boundary  of  the  first  knife  edge.  Accompanying  the  solid  line  is  a 
schematic  representation  of  a  curve  of  constant  attenuation  (shown  as  a  dashed  line) 
corresponding  to  the  same  excess  loss  as  that  evaluated  at  the  first  knife  edge.  Field 
points  falling  within  sector  (0)  are  assigned  a  partial  excess  loss  equal  to  the  standard 
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smooth  earth  attenuation  loss,  provided  it  does  not  exceed  the  value  of  the  upper  bound 
contour,  and  are  otherwise  assigned  the  value  equal  to  that  of  the  upper  bound.  Thus  field 
points  lying  above  the  upper  bound  contour  are  assigned  the  standard  smooth  earth 
excess  loss,  whereas  field  points  below  the  contour  are  assigned  a  constant  excess  loss 
which  is  the  value  defining  the  upper  bound  contour. 

As  the  field  point  approaches  the  first  sector,  (which  is  dominated  by  single  knife  edge 
diffraction  involving  the  transmitter  and  the  terrain  obstruction  labeled  KE1  in  the 
figure),  the  contribution  to  the  loss  from  the  first  knife  edge  is  assigned  a  value  of  zero 
until  the  clearance  ratio  above  the  first  knife  edge  exceeds  a  value  of -0.78.  The  geometry 
corresponding  to  this  transitional  value  of  the  clearance  ratio  is  indicated  in  the  figure  by 
a  square  bracket  labeled  Tl.  (In  addition,  the  geometry  corresponding  to  a  clearance  ratio 
of  -.78  is  indicated  for  each  sector  by  a  like-numbered  square  bracket  slightly  above  the 
respective  shadow  boundary.)  Field  points  below  the  critical  clearance  ratio  of  transition 
region  Tl  are  assigned  partial  losses  given  by  the  standard  Fresnel  integral  for  a  perfectly 
absorbing  knife  edge  [8],  subject  to  an  upper  bound  given  by  the  loss  computed  at  the 
next  higher  geodesic  point.  The  value  of  the  upper  bound  loss  for  each  sector  is  shown  in 
the  figure  as  a  dashed  contour  line  near  the  lower  bound  of  the  sector  and  intersecting  the 
corresponding  following  knife  edge.  As  the  vertical  descent  is  continued  a  non-zero  value 
of  excess  partial  loss  is  contributed  by  each  sector  encountered,  provided  the  clearance 
ratio  exceeds  -.78  and  provided  that  it  does  not  exceed  this  upper  attenuation  limit. 

As  the  vertical  descent  enters  the  lowermost  sector(4),  the  total  excess  path  loss  has 
accumulated  a  partial  excess  loss  from  each  of  the  four  sectors  given  by  the  upper 
attenuation  limit  of  that  sector  and  the  final  partial  loss  is  that  due  to  the  fourth  knife  edge 
with  effective  source  at  the  third  diffraction  center.  The  following  section  will  examine  a 
more  general  terrain  configuration  including  a  significant  smooth  earth  component  which 
will  require  elementary  smooth  earth  diffraction  methods  to  compute  the  corresponding 
partial  diffraction  contribution. 

3.2.2.  Discussion  of  irregular  terrain  including  portions  of  uniform 
elevation  (smooth  earth) 

A  more  general  terrain  profile  including  an  optically  significant  portion  of  uniform 
elevation  appears  in  Figure  2  and  is  further  characterized  by  two  non  inter- visible  knife 
edges.  The  diffraction  sectors  in  the  figure  include  that  marked  (0)  above  the  shadow 
boundary  of  the  first  knife  edge  (KE1),  sector  (1)  below  the  first  shadow  boundary  and 
above  the  first  tangent  line  passing  through  point  TNI,  sector  (2)  modeled  by  smooth 
earth  diffraction,  between  the  first  and  second  tangent  lines  passing  respectively  through 
points  TNI  and  TN2,  and  finally  that  marked  (3)  which  occupies  the  region  below  the 
shadow  boundary  of  the  second  knife  edge. 
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Figure  2.  General  Terrain  Profile  with  Two  Non-Intervisible  Knife  Edges 

An  objective  in  the  design  of  the  current  model  is  that  it  reduces  to  smooth  earth 
attenuation  in  the  absence  of  irregular  terrain  features.  In  order  to  accomplish  the 
reduction  to  smooth  earth  attenuation  in  the  absence  of  irregular  terrain,  field  points 
located  in  sector  zero  are  understood  by  default  to  sustain  an  initial  loss  due  only  to 
smooth  earth  attenuation.  As  the  transition  to  sector  one  is  approached,  the  smooth  earth 
attenuation  loss  is  restricted  to  an  upper  bound  given  by  its  value  at  the  first  diffraction 
center.  Continuity  across  the  sector  boundary  is  accomplished  by  defining  the  smooth 
earth  attenuation  as  its  conventional  functional  value  if  less  than  the  upper  bound  loss  and 
to  otherwise  defines  the  loss  as  equivalent  to  the  upper  bound  loss.  Consequently,  a 
partial  excess  loss  given  by  the  standard  functional  form  is  assigned  for  field  points  lying 
above  the  transition  contour  and  points  lying  below  the  transition  contour  are  assigned  a 
partial  loss  equivalent  to  the  its  defining  constant  value. 

As  the  boundary  between  sectors  0  and  1  is  approached,  an  additional  loss  due  to  the 
knife  edge  labeled  KE1  becomes  significant  as  the  clearance  ratio  above  the  first  knife 
edge  approaches  the  transition  value  of  -0.78.  The  geometry  along  the  second  vertical 
line  corresponding  to  this  transition  geometry  is  labeled  T1  in  Figure  2.  The  first  sector 
is  bounded  below  by  the  line  passing  through  point  TNI  which  is  the  point  of  tangency 
with  the  smooth  terrain.  The  contour  corresponding  to  the  limiting  value  of  loss 
contributed  by  the  first  knife  edge  is  schematically  shown  as  a  curved  dotted  line 
somewhat  above  the  tangency.  Below  the  transition  region  (T 1  )and  within  sector  one,  the 
partial  loss  contributed  by  the  first  knife  edge  is  given  by  the  conventional  functional 
form  if  less  than  this  upper  bound  value  which  is  also  the  limiting  value  assigned  for 
those  losses  otherwise  exceeding  the  limit. 


7 


Continuing  the  vertical  descent,  as  field  points  approach  the  tangent  line  through  TNI,  an 
additional  loss  attributable  to  attenuation  by  the  encroaching  smooth  earth  is  expressed 
for  points  below  the  smooth  earth  transition  region,  labeled  T2  in  Figure  2.  (The  choice 
to  use  smooth  earth  attenuation  is  justified  by  the  existence  of  smooth  terrain  between 
knife  edges  KE1  and  KE2  and  is  thus  modeled  with  an  effective  transmitter  at  the  same 
height  as  the  first  knife  edge  and  approaching  a  limiting  value  given  by  the  loss  evaluated 
at  the  second  knife  edge.)  The  additional  smooth  earth  loss  increases  in  value  in  the 
descent  through  sector  2  until  it  reaches  its  limiting  value  corresponding  to  the  dotted 
contour  above  the  shadow  boundary  of  the  second  knife  edge  (KE2).  Field  points  below 
the  upper  bound  contour  are  assigned  a  partial  diffraction  loss  equal  to  that  defining  the 
upper  bound.  A  third  partial  diffraction  loss  due  to  the  second  knife  edge  (KE2)  is 
sustained  by  those  field  points  lying  below  the  corresponding  transition  region  T3  in 
Figure  2.  This  loss  is  supplemented  by  an  additional  shadowing  loss  modeled  as  a  smooth 
earth  attenuation  of  a  previously  knife  edge  diffracted  field  as  discussed  in  the  appendix 
(13). 

3.3.  Extraction  of  Geodesic  Paths 

Since  the  diffraction  centers  essential  to  the  calculation  of  the  path  loss  are  confined  to 
local  geometrical  convexities  which  define  the  geodesic  path  between  relevant  antennas, 
it  becomes  necessary  to  extract  the  relevant  geodesic  geometry  intrinsic  to  the  terrain 
profile.  This  section  will  describe  a  diffraction  center  identification  process  intended  to 
avoid  a  possible  oversight  of  similar  methods  which  may  ignore  terrain  features  not 
among  the  vertices  of  the  dominant  geodesic  path  but  optically  relevant  to  the  calculation. 

Although  the  dominant  diffraction  centers  are  the  vertices  of  the  geodesic  path  between 
transmitter  and  receiver,  the  loss  calculation  may  be  significantly  affected  by  sub¬ 
geodesic  points,  especially  those  within  a  Fresnel  zone  of  the  dominant  geodesic  ray 
paths.  Restricting  the  optically  relevant  environment  to  the  external  geodesic  path  may 
lead  to  discontinuities  in  path  loss  calculation  because  the  evolution  of  the  terrain  may 
introduce  dramatic  alterations  of  the  points  comprising  the  dominant  geodesic  path, 
permitting  an  erroneous  neglect  of  optically  significant  terrain  features  no  longer  among 
the  dominant  geodesic  vertices.  Inclusion  of  sub-geodesic  paths  assures  that  points 
originally  dominant  in  an  evolving  terrain  profile  may  continue  to  exert  their  influence  as 
the  path  is  extended. 


As  an  example,  consider  the  terrain  profile  appearing  in  Figure  3  which  consists  of  a 
random  assortment  of  five  knife  edges  distributed  over  a  25  km  path.  Critical  terrain 
points  of  this  profile  have  been  labeled  by  the  order  they  appear  within  the  profile. 
Terrain  points  labeled  (0),  (11),  and  (22),  marked  with  red  arcs,  are  the  dominant 
geodesic  points  for  the  first  10  km  of  the  path,  but  point  (22)  is  excluded  from  the  new  set 
of  dominant  geodesic  points  which  results  if  the  path  is  extended  to  point  (27)  which 
replaces  point  (22)  as  the  third  point  of  the  resulting  dominant  path  marked  with  green 
circles. 

A  branching  of  geodesic  paths  is  implied  from  point  1 1  of  the  first  geodesic  path  to  point 
27  which  is  the  first  unique  point  of  the  second  geodesic  path,  marked  in  green,  which  is 
then  extended  to  point  33.  This  structure  of  the  dominant  geodesic  path  remains  intact  for 
paths  of  a  lesser  range  coordinate  than  profile  point  38  where  a  second  branching  from 
point  27  occurs  to  access  the  first  unique  point  of  the  outermost  geodesic  path  marked 
with  blue  arcs.  Since  the  first  set  of  optical  paths  marked  with  red  and  green  arcs  conform 
to  the  definition  of  optically  relevant  points  for  the  shorter  paths,  a  generalization  of 
optical  relevance  would  include  all  points  considered  optically  relevant  anywhere  along 
an  exhaustive  traversal  of  the  entire  path  profile.  Defining  points  contained  in  a 
continuous  geodesic  (convex)  path  originating  at  the  transmitter  as  elements  of  a  geodesic 
seam  regardless  of  length,  the  optically  relevant  points  may  be  defined  as  all  points  lying 
on  the  geodesic  seams  generated  by  the  specified  invariant  transmitter.  Thus  each  of  the 
three  paths  discussed  above  originate  at  the  transmitter,  define  a  locus  of  continuous 
convexity,  comprise  a  geodesic  seam  and  are  therefore  included  among  a  universal  set  of 
optically  relevant  points,  as  defined  here. 

The  geodesic  seams  are  located  by  recording  the  indices  of  the  geodesic  path  to  a 
specified  point  on  the  path  profile  as  a  function  of  the  range  coordinate  of  that  point.  The 
geodesic  paths  are  characterized  by  monotonically  decreasing  slopes  and  are  therefore  of 
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a  negative  second  derivative.  The  relevant  geometry  for  geodesic  path  extension 
inspection  is  depicted  in  Figure  4,  where  it  is  understood  that  final  two  points  of  a  known 
geodesic  path,  labeled  10  and  1 1 ,  have  already  been  established.  Inspection  of  the  point  at 
12  for  its  suitability  as  an  extension  of  the  known  geodesic  path  also  requires  knowledge 
of  the  point  just  beyond  point  12  and  shown  as  13  in  the  figure.  The  continuity  of 
negativity  requires  that  the  curvatures  of  the  approach  and  exit  vectors  associated  with 
point  12  are  both  negative.  (The  vector  r01  from  point  (0)  terminating  at  point  (1)  can  be 

used  to  define  mutually  orthogonal  unit  vectors  e\  =ro\ and  ej  as  kxe j .  The  orientation 
of  r\ 2  with  respect  to  fgi  is  negative  in  the  formal  sense  that  r\2'e2  <  ^  •  Alternatively 

the  sign  of  the  curvature  may  be  defined  as  the  sin  of  angle  Z(I0-I1-I2)  which  is  critically 
negative  if  it  exceeds  n )  This  requires  that  the  approach  angle  to  (12),  Z(I0-I1-I2)  , 
terminating  on  the  point  in  question,  and  the  exit  angle  Z(I1 -12-13),  centered  on  the  point 
in  question,  must  both  be  negative.  Confirmation  of  both  negative  slopes  justifies  the 
inclusion  of  the  point  (12)  as  a  convex  extension  of  the  geodesic  path  whose  current 
leading  element  is  point  II. 


For  computational  purposes,  the  proper  curvatures  of  approach  and  exit  are  confirmed  for 
each  additional  extension  point  and  the  process  is  repeated  until  a  point  is  found  which 
violates  the  continuity  of  negativity.  If  the  point  in  question  is  found  to  exhibit  a  negative 
approach  angle  but  a  positive  exit  angle,  it  is  recognized  as  being  the  non  convex 
terminus  of  the  current  geodesic  path  and  is  discarded  from  further  consideration.  If  the 
point  however  is  found  to  exhibit  a  positive  approach  angle,  the  current  path  is  terminated 
but  the  inspected  point  is  retained  for  possible  insertion  into  a  higher  order  geodesic  path. 
The  laminated  structure  of  the  geodesic  seams  assures  that  a  new  geodesic  path  to  such 
an  excluded  point  may  be  established  from  a  point  within  the  current  outermost  geodesic 
path. 
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As  exemplified  in  Figure  3,  the  first  (red)  geodesic  path  terminates  at  point  (22)  and  a 
new  geodesic  path  is  initiated  at  point  (27)  as  an  extension  from  previous  point  (11) 
which  is  the  last  point  among  those  of  the  previous  geodesic  path  to  satisfy  the  double 
negative  curvature  constraint  for  points  (0,  11  27,  28)  serving  as  (10,  II,  12,  13) 
respectively.  Points  1-10,  12-21  and  23-26  are  eliminated  in  this  process  because  none 
are  properly  configured  to  provide  negative  exit  curvature  when  considered  as  extensions 
of  any  previous  geodesic  point.  Continuing  the  analysis  beyond  point  27,  points  28-32  are 
eliminated  because  , although  each  as  point  12  admits  to  a  negative  approach  curvature 
from  points  (10,  II;  12,  12  +  1)  =  (11,  27;  12,  12  +  1),  none  is  properly  configured  with 
respect  to  the  corresponding  following  point  at  13  =  12  +  1  to  provide  the  requisite 
negative  exit  curvature.  An  extension  is  found  at  point  33  which  exhibits  the  proper 
approach  and  exit  curvatures  as  defined  by  the  vector  relations  of  point  (10,  II;  12, 13)  = 
(11,  27;  33,  34)  respectively.  No  suitable  extensions  among  the  points  34-37  are  found 
because,  although  each  exhibits  the  desired  approach  curvature,  each  also  exhibits  an 
undesirable  exit  curvature  when  approached  form  point  33.  Since  the  approach  curvature 
to  point  38  from  {10,  II;  12  13}  =  (27,  33;  38,  39}  is  positive,  the  geodesic  path  is 
terminated,  a  new  branch  point  from  the  previous  geodesic  path  at  point  27  is  identified 
and  the  third  and  final  geodesic  path  is  extended  beyond  point  38  to  its  final  extension  at 
point  46. 

The  methods  described  above  for  irregular  terrain  are  directly  applicable  to  portions  of 
the  profile  determined  to  be  of  uniform  elevation  as  expected  at  sea  level.  The  upper  right 
panel  of  Figure  5  exemplifies  a  terrain  profile  consisting  of  irregular  features,  manifested 
here  as  knife  edges  at  points  (318)  and  (399)  and  of  an  extended  portion  of  uniform 
elevation  between  the  transmitter  at  point  (0)  and  point  (315).  Extension  of  the  geodesic 
path  beyond  point  0  involves  examining  approach  and  exit  slopes  of  (10,  II,  12, 13)  for  12 
beyond  0.  (Points  10  and  II  are  formally  initialized  as  lying  at  the  base  and  apex  of  the 
transmitter  which  is  assigned  terrain  index  0)  Points  1-108  are  excluded  because  they 
each  exhibit  positive  exit  slope  despite  the  negative  approach  curvature  from  point  0. 
Point  109  is  the  first  point  beyond  the  smooth  earth  horizon  of  point  (0)  and  is  also  the 
first  point  to  exhibit  a  negative  exit  curvature  in  the  approach  from  point  0  and  therefore 
the  first  extension  of  the  geodesic  path.  The  ray  path  from  point  (0)  to  its  horizon  at  point 
(109)  may  be  described  as  the  tangential  insertion  of  the  ray  path  from  the  transmitter  at 
point  (0)  into  the  geodesic  path  extended  by  a  smooth  earth  arc  beginning  at  point  109. 
When  restricted  to  uniformly  elevated  terrain  beyond  an  elevated  point  of  a  ray  path,  the 
search  for  path  extensions  with  continuous  negative  curvature  has  then  led  to  the 
detection  of  the  horizon  of  the  elevated  terrain  point. 
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Figure  5.  Smooth  Earth  with  Inter-Visible  Knife  Edge 


Upon  tangential  insertion  into  the  uniformly  elevated  portion  of  the  profile,  here  at  the 
horizon  point  (109),  all  subsequent  points  satisfy  the  criterion  of  extended  geodesic 
curvature  until  irregular  terrain  is  encountered  at  point  (318).  Responding  as  described 
earlier  to  a  detection  of  obstructed  geodesic  geometry,  a  search  is  conducted  among  the 
established  geodesic  points  for  a  branch  point  to  this  first  unique  point  of  the  next  higher 
order  geodesic  path.  From  an  inspection  of  points  previous  to  point  318  for  the  proper 
approach  and  exit  slopes,  it  becomes  apparent  that  point  0  serves  as  the  branch  point  to 
the  next  higher  order  geodesic  path  and  that  point  318  is  the  first  unique  point  of  the  path 
which  has  been  marked  in  green  in  the  figure.  The  branching  point  is  at  point  (0)  because 
the  knife  edge  at  point  (318)  is  above  the  projected  horizon  of  point  0  shown  in  the 
Figure  5  as  an  extended  tangent  intersecting  the  transmitter  horizon  at  point  (109). 

A  variation  of  this  geometry  is  shown  in  Figure  6  in  which  the  first  knife  edge  is  now 
below  the  projected  horizon  of  the  transmitter.  As  with  inter-visible  obstructions 
separated  by  smooth  earth  discussed  earlier,  the  first  geodesic  path  proceeds  from  the 
transmitter  at  point  (0)  to  its  horizon  which  for  this  profile  is  now  at  point  127  and  along 
the  smooth  earth  arc  to  the  first  obstruction  at  point  396.  The  first  geodesic  path  is 
terminated  and  a  search  is  conducted  for  proper  branching  to  a  higher  order  geodesic  path 
as  before.  It  becomes  apparent  from  the  resulting  inspection  of  points  that  point  (219), 
beyond  the  first  tangential  insertion  at  point  (127)  and  also  a  point  within  the  uniformly 
elevated  portion  of  the  terrain,  is  the  last  point  among  those  inspected  which  satisfies  the 
curvature  requirements  for  geodesic  extension.  The  branching  from  point  219  is  the 
greater  terminus  of  the  smooth  geodesic  arc  between  points  127  and  219,  each  of  which  is 
a  tangential  insertion  focus  and  local  (smooth  earth)  horizon  for  the  corresponding 
irregular  terrain  features. 
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Figure  6.  Knife  Edge  Beyond  Horizon 


3.4.  Designation  of  Diffraction  Centers  and  Forward 
Propagation  Shadows 

Provided  the  geodesic  seams  of  the  terrain  profile,  the  diffraction  centers  may  be 
identified  according  to  a  methodology  described  in  this  section.  It  is  customary  to 
minimally  require  that  terrain  vertices  to  be  modeled  as  diffraction  centers  must  block  the 
direct  Line  of  sight  ray  path  between  an  emitter  and  other  portions  of  the  terrain.  The 
geodesic  seams  identified  in  the  previous  section  consist  of  serialized  lines  segments 
between  adjacent  vertices  of  the  optical  ray  path  and,  for  those  paths  encountering 
smooth  earth  profiles,  continuous  geodesic  arcs  of  uniform  elevation  bounded  by 
transitional  tangent  rays  from  both  endpoints  to  the  corresponding  irregular  terrain 
relation. 

The  geodesic  seams  provide  the  vertices  of  the  discrete  portions  of  the  optical  ray  as  well 
as  any  continuous  geodesic  segments  introduced  by  portions  of  terrain  determined  to  be 
of  uniform  elevation.  Since  the  points  of  the  geodesic  path  necessarily  exhibit  continuity 
of  negative  curvature,  they  also  satisfy  the  earlier  requirement  of  successive  ray  path 
obstruction.  The  vertices  of  the  discrete  portion  of  the  geodesic  seams  appear  to  then 
qualify  as  centers  of  diffraction  which  for  the  purposes  of  this  report  are  modeled  as 
perfectly  absorbing  knife  edges.  The  continuous  portions  of  the  ray  paths,  as  exemplified 
by  the  arc  between  points  (127)  and  (219)  in  Figure  6,  and  their  most  immediate  irregular 
relations,  as  exemplified  by  points  (0)  and  (396)  in  the  Figure  6,  define  elements 
consisting  of  effective  antennas  of  finite  height  separated  by  uniformly  elevated  terrain 
and  may  therefore  be  modeled  as  localized  smooth  earth  attenuation  partial  diffraction 
factors. 

According  to  the  identification  procedure  thus  described,  the  diffraction  centers  of  the 
example  profiles  shown  in  Figure  7,  Figure  9,  and  Figure  1 1  are  identified  respectively  in 
Figure  13  -  Figure  15  in  which  the  knife  edge  diffraction  centers  of  the  discrete  geodesic 
portion  are  indicated  by  the  knife  edge  attenuation  label  AKE(IL,  10,  IG)  where  the  three 
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integer  arguments  reference  the  ordering  of  points  in  the  geodesic  path  serving  as 
effective  source,  diffraction  center  and  limiting  field  point,  respectively.  The  smooth 
earth  attenuation  diffraction  centers  are  indicated  in  Figure  14  and  Figure  15  by  the 
smooth  earth  attenuation  label  ASE(I0,  II,  12, 13)  in  which  the  integer  arguments  indicate 
geodesic  indices  of  points  serving  respectively  as  source,  first  and  second  tangents,  and 
limiting  attenuation  field  point.  The  effective  transmitter  geometry  for  the  partial 
diffraction  losses  due  to  the  knife  edge  at  point  (II-3)  is  defined  by  the  geodesic  path  just 
preceding  the  knife  edge,  which  in  Figure  15  is  point  (II-2).  Partial  losses  contributed  by 
this  knife  edge  are  then  computed  with  an  effective  transmitter  located  at  point  (II-2) 
which  serves  as  the  tangential  origin  of  the  ray  path  from  the  smooth  earth  to  the  knife 
edge  at  point  (II-3). 
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Figure  8.  TMEM  Propagation  Calculation  at  20.08  MHz 
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Figure  10.  TMEM  Smooth  Earth  Inter-Visible  Knife  Edge  at  20.08  MHz 
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Figure  12.  TMEM  Knife  Edge  Beyond  Horizon  at  20.08  MHz 


RANGE  COORDINATE  (KM) 

Figure  13.  Diffraction  Centers  for  Sample  Terrain  Profile 


RANGE  COORDINATE  (KM) 

Figure  14.  Smooth  Earth  Tangency  with  Visible  Discrete  Diffraction  Center 
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Figure  15.  Knife  Edge  Beyond  Smooth  Earth  Horizon 


3.5.  Ray  Path  Construction  and  Total  Loss  Calculation  for 
Arbitrary  Field  Point 

With  the  geodesic  seams  of  the  terrain  profile  established  as  the  essential  intermediates 
for  ray  path  construction,  it  becomes  possible  to  establish  the  geodesic  path  and  therefore 
the  discrete  and  continuous  diffraction  centers  of  relevance  to  an  arbitrary  field  point.  As 
in  the  analysis  associated  with  geodesic  path  extension,  the  geodesic  path  to  an  arbitrary 
field  point  may  be  regarded  as  the  composite  of  the  geodesic  path  to  the  critically  convex 
terrain  point,  defined  below,  and  the  ray  path  from  that  point  to  the  field  point.  Referring 
again  to  the  notation  and  geometry  of  Figure  16,  the  convex  critical  point  is  the  terrain 
point  (II)  of  largest  range  with  most  immediate  lesser  relation  (10)  of  the  geodesic  path 
which  admits  to  a  negative  curvature  in  the  approach  to  the  field  point  (12).  Associated 
with  the  convex  critical  vertex  is  a  radial  sector  bounded  above  by  the  extension  of  the 
geodesic  path  approaching  the  vertex  and  bounded  below  by  extension  of  the  path  to  the 
next  geodesic  point  of  the  path.  As  shown  in  Figure  18  the  sectors  may  be  further 
subdivided  by  intermediate  obstructions  if  the  defining  knife  edge  also  serves  as  a  branch 
to  lesser  order  geodesic  paths  interior  to  the  outer-most  geodesic  path.  In  the  following 
section  the  more  general  aspects  of  shadow  boundary  continuity  associated  with  traversal 
of  various  diffraction  sectors  including  sub  geodesic  obstructions  will  be  examined  from 
a  more  comprehensive  perspective. 
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3.5.1.  Continuity  of  attenuation  across  shadow  boundaries 

Provided  the  geometry  of  diffraction  centers,  associated  upper  bound  shadow  boundaries 
and  range  dependent  lower  bound  restriction  of  any  sub  geodesic  obstructions,  this 
section  will  examine  the  rationale  for  enforcement  of  total  attenuation  continuity  as  the 
diffraction  sectors  are  traversed  and  the  total  loss  is  accumulated.  Before  considering  the 
more  general  terrain  geometry  in  which  a  dominant  geodesic  path  is  obstructed  by  a 
number  of  sub  geodesic  features,  the  basics  of  shadow  boundary  continuity  enforcement 
may  be  described  as  applied  for  the  simple  truncated  wedge  appearing  in  Figure  17.  Each 
of  the  two  diffraction  centers  which  include  the  vertex  at  the  apex  of  the  wedge  and  the 
vertex  to  the  right  of  the  apex  is  characterized  by  one  of  two  mutually  adjacent  diffraction 
sectors.  The  upper  bound  of  the  apical  diffraction  sector  is  its  shadow  boundary  which 
appears  in  the  figure  as  a  ray  path  extended  beyond  the  vertex  with  the  same  slope  as  the 
approach  to  the  vertex  from  the  transmitter.  The  lower  bound  of  the  apical  diffraction 
sector  is  the  radial  extension  of  the  path  from  the  vertex  to  the  next  diffraction  center.  As 
described  below,  a  diffraction  ray  path  and  the  associated  partial  diffraction  loss  will  be 
computed  for  any  unobstructed  field  point  within  the  respective  diffraction  sector. 
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RANGE  COORDINATE  (KM) 

Figure  17.  Shadow  Boundary  Continuity  for  Truncated  Wedge 


Along  the  vertical  line  passing  through  point  A,  the  smooth  earth  attenuation  is  evaluated 
from  the  upper  margin  to  the  intersection  with  the  dashed  curved  contour  passing  through 
point  A  corresponding  to  its  value  at  the  first  knife  edge  and  is  assigned  this  constant 
value  below  this  contour.(The  assignment  of  the  upper  attenuation  limit  above  the 
shadow  boundary  assures  continuity  with  the  loss  as  computed  below  the  shadow 
boundary  as  modeled  by  Epstein  and  Peterson  [1]  who  express  the  total  loss  as  the  sum  of 
two  independent  single  knife  edge  losses  including  that  due  to  the  first  knife  edge  with 
the  second  knife  edge  as  its  effective  field  point  and  that  due  to  the  second  knife  edge 
with  the  first  knife  edge  serving  as  its  effective  source.) 

Along  the  vertical  line  containing  points  B  and  C,  the  first  knife  edge  attenuation  is 
assigned  the  value  of  1  (zero  loss)  above  transition  region  B  where  the  scaled  clearance 
argument  reaches  a  value  of  -.78  The  knife  edge  attenuation  is  evaluated  by  the 
associated  knife  edge  attenuation  function  based  on  a  ID  Fresnel  integral  between 
transition  region  B  and  point  C  where  is  reaches  its  maximum  attenuation  value 
corresponding  to  the  dashed  curved  contour  and  is  assigned  this  constant  value  below 
point  C.  Along  the  vertical  line  passing  through  transition  region  D,  the  second  knife 
edge  attenuation  is  assigned  a  value  of  1  above  transition  region  D  where  the  scaled 
clearance  argument  attains  a  value  of  -.78  and  is  assigned  the  value  of  the  associated 
knife  attenuation  function  based  on  a  ID  Fresnel  integral  value  below  transition  region  D. 
A  more  general  diffraction  center  obstructed  by  sub-geodesic  terrain  obstacles  is 
considered  in  the  following  section. 
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3.5.2.  Discussion  of  Generalities  Introduced  by  Sub  Geodesic 
Obstructions  (profiles  with  shadows  obstructions) 

The  more  general  terrain  profile  including  sub  geodesic  obstructions  of  a  lower  order 
than  the  external  geodesic  path  is  shown  in  Figure  18.  Apparent  in  the  figure  are  the 
complications  introduced  by  optical  obstructions  of  a  diffracted  ray  path  and  the 
geometrical  considerations  for  inclusion  of  a  diffraction  component  from  a  diffraction 
center  not  among  those  of  the  dominant  geodesic  path.  In  particular,  the  diffracted  rays 
attributable  to  the  second  knife  edge  (KE2)  as  illuminated  by  the  transmitter,  point  (0), 
encounter  range-dependent  obstructions  (KE(3-4))  between  the  second  knife  edge  and  the 
dominant  geodesic  outlet  at  KE5.  The  existence  of  an  unobstructed  ray  path  to  a  field 
point  between  KE2  and  KE5  can  be  confirmed  only  by  establishing,  for  each  of  the  three 
subintervals  imposed  by  the  intervening  knife  edges  that  the  field  point  in  question  lies 
above  the  local  shadow  boundary  shown  as  a  solid  line  in  each  subinterval.  Beyond  KE5 
the  knife  edge  labeled  KE2  is  no  longer  among  the  points  in  the  dominant  geodesic  path 
extending  directly  from  the  transmitter  to  KE5.  The  partial  diffraction  loss  due  to  this 
secondary  diffraction  center  may  be  described  as  a  function  of  elevation  in  a  descent 
along  the  vertical  line  appearing  in  the  figure  between  knife  edge  KE5  and  the  right 
margin.  Between  the  upper  margin  and  transition  region  T  (shown  as  a  square  bracket 
above  the  shadow  boundary),  the  partial  diffraction  loss  due  to  (KE2)is  the  standard 
single  knife  edge  loss  involving  source  at  point  (0),  diffraction  center  at  point  (KE2)  and 
limiting  field  point  at  point  (KE5)  for  which  the  notation  AKE  (025)  is  used  in  the  figure. 
The  limiting  attenuation  due  to  this  diffraction  center  is  defined  by  its  value  at  KE5. 
Below  the  transition  region,  the  partial  diffraction  loss  due  to  the  second  knife  edge 
beyond  KE5  is  assigned  the  constant  value  given  by  this  upper  bound  attenuation.  The 
similarity  of  this  modeling  to  that  used  for  diffraction  centers  within  the  dominant 
geodesic  path  assures  a  continuous  transition  to  geometries  in  which  the  considered 
obstruction  would  be  among  the  dominant  diffraction  centers. 


RANGE  COORDINATE  (KM) 

Figure  18.  Sub  Geodesic  Diffraction 
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3.6.  Fused  wedges 

In  anticipation  of  the  comprehensive  desirability  of  a  propagation  model  which  is  at  once 
capable  of  treating  terrain  including  both  irregular  and  artificial  features  with  piecewise 
linear  faces  as  in  polygonal  obstacles,  it  becomes  necessary  to  include  a  terrain  geometry 
recognition  feature  which  detects  portions  of  piecewise  linear  terrain  including  those 
features  essential  in  the  processing  of  the  diffraction  centers.  A  representative  profile 
including  both  irregular  terrain  and  a  block  shaped  obstacle  intended  to  suggest  a 
schematic  representation  of  a  building  appears  in  Figure  19. 


Figure  19.  Terrain  Profile  with  Polygonalizable  and  Irregular  Terrain 

Since  the  defining  feature  of  the  artificial  obstacles  is  the  presence  of  any  piecewise 
geometrically  linear  faces  which  can  only  be  represented  with  discrete  profile 
assignments  by  collinear  points,  the  algorithm  described  earlier  for  identifying  geodesic 
seams  imbedded  in  the  terrain  profile  was  enhanced  to  include  a  local  co-linear  detection 
feature  which  in  its  examination  of  terrain  points  simultaneously  inspects  for  co-linearity 
involving  each  point  in  the  profile  and  the  corresponding  two  previous  contiguous  points. 
Upon  detecting  any  isolated  three  such  co-linear  points,  a  line  segment  is  initiated  and 
possibly  extended  by  any  subsequent  points  conforming  to  the  co-linearity  criterion.  The 
relevant  line  segment  is  extended  until  a  non  linear  terminus  is  detected  and  the 
continuous  set  of  points  comprising  a  co-linear  segment  is  provided  an  identifying 
segment  index  and  the  constituent  points  are  each  labeled  as  components  of  the  segment. 
Anticipating  the  use  of  single  and  multiple  wedge  diffraction  according  to  the  methods  of 
the  uniform  geometrical  theory  of  diffraction  (UTD)  to  model  the  diffraction  associated 
with  the  vertices  defined  by  adjacent  co-linear  segments,  the  dependence  of  UTD  upon 
only  the  internal  wedge  angle  defined  by  adjacent  co-linear  segments  permits  the 
elimination  of  the  interstitial  points  of  the  detected  line  segments  in  which  only  the 
endpoints  are  retained.  Since  any  two  adjacent  co-linear  segments  will  have  a  common 
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point,  the  composite  construction  consisting  of  the  two  adjacent  line  segments  and  their 
common  point  vertex  is  sufficient  to  define  a  UTD  diffraction  center  provided  the 
geodesic  path  to  the  diffraction  center  is  also  provided.  As  an  example,  the  profile 
appearing  in  Figure  19  is  observed  to  contain  a  continuous  sequence  of  co-linear  points 
defining  the  upper  face  of  a  rectangular  obstacle  and  elimination  of  the  associated 
interstitial  points  produces  the  (co-linearly)  reduced  profile  appearing  in  Figure  20. 


Figure  20.  Terrain  Profile  with  Polygonalizable  and  Irregular  Terrain  After  Removal 

of  Interstitial  colinearities  (Reduced  Terrain) 

Since  the  endpoints  of  the  artificial  terrain  feature  are  also  diffraction  centers,  the  two 
complementary  faces  for  completing  the  geometry  of  the  two  (fused)  wedges  are  both 
taken  as  their  most  immediate  non  collinear  relations.  Thus  the  co-linearly  reduced 
profile  appearing  in  Figure  20  includes  the  two  diffraction  centers  labeled  in  the  figure 
and  the  critical  internal  wedge  angles  required  for  wedge  diffraction  are  defined  as  the 
angles  spanned  in  a  rotation  about  the  diffraction  centers  between  the  corresponding 
bounding  collinear  faces.  In  Figure  20  the  two  diffraction  centers  labeled  Qoi2  an<^  0123 

are  then  characterized  by  the  internal  wedge  angle  complementary  to  the  arc  centered  at 
the  labeled  point  and  extending  between  its  most  immediate  (reduced)  terrain  relations. 
The  identification  of  diffraction  centers  defined  by  two  co-linear  segments  is 
accomplished  by  assigning  a  value  of  one  to  a  UTD  identification  flag  which  is  assigned 
for  all  diffraction  centers  and  otherwise  initialized  as  zero.  Any  subsequent  calculation  of 
diffraction  related  quantities  involving  such  diffraction  centers  will  be  executed 
according  to  the  relevant  UTD  formalism  including  single  and  multiple  serialized  and 
fused  wedge  composite  structures.  The  following  section  will  describe  utilization  of 
knowledge  of  the  geodesic  paths  to  alls  diffraction  centers  and  their  forward  propagation 
direction  optical  obstructions  to  extract  the  relevant  arguments  for  UTD  diffraction. 
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Since  the  geodesic  paths  and  forward  propagation  optical  shadows  are  available  for  all 
diffraction  centers  including  those  with  positive  UTD  identification  flags  and  since 
geodesic  illumination  is  assumed  for  all  relevant  diffraction  calculations,  the  standard 
source  coordinates  p'  and  <j>’  are  the  polar  coordinates  of  the  source  point  associated 
with  the  most  immediate  geodesic  relation  in  the  direction  of  the  transmitter.  These 
source  coordinates  are  computed  in  a  coordinate  system  with  the  diffraction  center  as 
origin  and  with  rotated  Cartesian  axes  defined  by  the  unit  vector  el  in  the  direction  of  the 
most  immediate  terrain  relation  in  the  direction  of  the  source  and  therefore  aligned  along 
the  illuminated  side  of  the  wedge,  and  an  associated  orthogonal  unit  vector  defined  by 

(e,  )x(&).  The  field  point  coordinates,  p  and  <f>,  are  the  radial  coordinates  of  the  relevant 
receiver  target  as  computed  in  the  same  coordinate  system.  The  processing  of  the 
diffraction  centers  designated  as  compatible  with  the  methods  of  UTD  is  expected  to 
proceed  identically  to  those  modeled  as  elements  of  serialized  knife  edge  or  smooth  earth 
diffraction  centers  including  establishment  of  immediate  local  source  relations  and 
evaluation  of  upper  bound  attenuations  at  any  optical  obstructions  in  the  direction  of 
forward  propagation.  The  specific  complex  partial  UTD  field  component  contributed  to  a 
sequence  of  serialized  diffractions  is  given  by  the  relation 
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where  E1q  1  is  the  complex  incident  field  as  propagated  along  the  relevant  geodesic  path 

to  the  diffraction  center  in  question,  D  is  the  UTD  diffraction  coefficient  of  Kouyoumjian 
and  Pathak  and  where  the  final  factor  is  the  spreading  factor  constructed  from  quantities 
st  and  s  which  are,  respectively,  the  sum  of  all  finite  ray  path  elements  encountered  in 
the  geodesic  ray  path  to  the  diffraction  center,  and  the  final  ray  path  component  from  the 
diffraction  center  to  the  field  point. 

In  the  interest  of  enhanced  conformity  with  earlier  methods  involving  diffraction 
simulation  with  perfectly  absorbing  knife  edges,  it  has  been  found  possible  to  use  this 
defining  basic  element  in  conjunction  with  image  theory  to  closely  approximate  the  loss 
patterns  of  wedge  diffraction.  Considering  the  evaluation  of  the  field  point  labeled  R  in 
Figure  21,  the  knife  edge  attenuated  fields  are  independently  evaluated  at  the  field  point 

r ,  to  obtain  EKE  (r) ,  and  at  the  virtual  image  r’  to  obtain  Eke (r1) .  The  virtual  receiver 
vector  r'  is  obtained  as  the  reflected  image  of  r  through  the  dark  surface  of  the  wedge 
(with  normal  unit  vector  e2  as  shown  in  Figure  21  according  to  the  reflection  operation: 

r'  =  r-2(r-e2]e2 
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Figure  21.  Use  of  Image  Theory  to  Simulate  Wedge  Diffraction  with  Fresnel  (Half 

Plane)  Diffraction  Theory 

A  reflection  coefficient  for  this  process  is  then  evaluated  consistent  with  a  field  points 
inclined  relative  to  the  reflecting  surface  with  skew  angle  y/  =  wx  -  <j) .  A  similar  pair  of 

partial  fields  is  then  constructed  as  before  using  the  image  of  the  transmitter,  obtained  as 
a  reflection  through  the  lit  face  of  the  wedge  which  also  requires  evaluation  of  the 
relevant  reflection  coefficient  with  an  inclination  from  the  reflecting  surface  of  skew 
angle  y/  =  <j> .  The  four  conjugate  partial  fields  are  then  superposed  to  obtain  the  total 
diffracted  field  as 


E  =  EKE  (rTX  ,rRx)+  RnEKE  (rTX  >rkx)  +  Ro[eKE  (rTX  >rRx)+  Rn^KE  {rTX  >rRX  )) 

The  agreement  with  the  exact  UTD  results  as  shown  in  the  vertical  excursion  plots  of 
Figure  22  and  the  contour  plots  of  Figure  23  and  Figure  24  is  encouraging. 
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Figure  23.  Propagation  Plane  Path  Loss  Contours  Obtained  Using  Knife  Edge 

Simulated  Wedge  Diffraction 


Figure  24.  Propagation  Plane  Path  Loss  Contours  for  Wedge  Diffraction  Obtained 

Using  Exact  (U)GTD 


4.  Results 

The  latest  trial  test  for  the  TMEM  propagation  algorithm  show  promising  results.  This 
section  compares  TMEM  to  other  existing  propagation  models  in  several  different  test 
cases. 

The  first  test  case  seeks  to  isolate  diffraction  from  the  other  phenomena.  In  this  case  a 
transmitter  at  1,278  meters  MSL,  vertical  polarization,  min  range  =  25  km,  max  range  = 
50  km,  min  height  =  1,248  meters,  max  height  =  3,248  meters  for  AGL  and  MSL, 
number  of  ranges  =  201,  number  of  heights  =  101,  the  terrain  profile  in  range  and  height 
in  meters  to  (0,  0),  (25, 1,  748),  and  (50, 0). 

The  PE  methods  demonstrate  the  underlying  physics  very  well.  TMEM  and  TIREM 
match  the  PE  methods  quite  well  in  the  shadow  region. 
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Figure  25.  Wedge,  Model  Comparisons  at  10  GHz,  Vertical  Polarization 


The  next  test  case  involves  an  80  km  terrain  segment  of  North  Table  Mountain,  obtained 
from  the  Institute  for  Telecommunication  Science.  In  this  case  a  transmitter  at  5  meters 
AGL,  vertical  polarization,  min  range  =  0  km,  max  range  =  80  km,  min  height  =  1,450 
meters,  max  height  =  2,200  meters  for  AGL  and  MSL,  number  of  ranges  =  301,  number 
of  heights  =  401. 

The  PE  method  demonstrates  the  underlying  physics  very  well.  TIREM  and  TMEM 
match  the  PE  methods  quite  well  in  the  shadow  region. 
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Figure  26.  ITS  Terrain  Data,  Model  Comparisons  at  10  GHz,  Vertical  Polarization 

The  last  test  case  involves  a  50  km  terrain  segment  of  the  Colorado  Mountains.  In  this 
case  a  transmitter  at  7  meters  AGL,  vertical  polarization,  min  range  —  0  km,  max  range  = 
50  km,  min  height  =  1,600  meters,  max  height  =  3,500  meters  for  AGL  and  MSL, 
number  of  ranges  =  251,  number  of  heights  =  301. 

The  PE  method  demonstrates  the  underlying  physics  very  well.  TMEM  and  TjpjgM 
match  the  PE  method  quite  well  in  the  shadow  region. 
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Figure  27.  Colorado  Mountains,  50  km  Range, 
Model  Comparisons  at  10  GHz,  Vertical  Polarization 
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5.  Further  Research 

Additional  research  planned  for  the  Terrain  Multipath  Electromagnetic  Model  is  detailed 
in  the  following  sections. 

5.1.  Adjustment  of  Excessive  Loss  for  Multiple  Knife  Edge 
Attenuations  and  for  Overlapping  Shadow  Boundaries 

The  essential  equation  for  computing  attenuation  as  modeled  by  a  perfectly  absorbing 
knife  edge  is  derived  from  asymptotic  properties  of  the  more  rigorous  formulations  which 
result  from  the  assumption  that  displacements  in  range  are  much  greater  than  the 
orthogonal  displacements  in  elevation.  The  exaggerated  excess  losses  observed  earlier  for 
the  proposed  model  in  comparison  with  TIREM,  displayed  again  in  Figure  28,  are  an 
illustration  of  a  violation  of  those  basic  assumptions  in  that  the  separation  of  knife  edges 
is  only  200  meters.  Although  somewhat  arbitrary  lower  limits  could  be  imposed  on  knife 
edge  spacing  ,  its  is  the  intention  of  this  effort  to  include  attenuation  compensation 
factors  which  will  properly  adjust  the  excessive  multiple  diffraction  resulting  from 
diffraction  centers  of  less  than  optimum  spacing.  Deygout  [88]  has  derived  an  expression 
which  serves  as  a  double  knife  edge  attenuation  correction  which  could  be  applied  for 
each  of  the  successive  knife  edges  pairings  defining  a  multiple  knife  edge  configuration 
and  is  of  the  specific  form 
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l  X) 
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where  q  and  p  are  scaled  knife  edge  heights  and  a  is  the  geometrical  adjacency 
parameter  derived  in  the  exact  two  knife  edge  study  of  Millington  et.  al.  [9].  If  the 
spacing  parameter  is  ideal,  a  assumes  the  value  of  k  /  2 ,  and  the  knife  edge  correction 
factor  is  1  but  as  alpha  approaches  a  vanishing  value,  the  correction  attenuation  factor 
approaches  2  for  a  loss  correction  that  results  in  6dB  less  loss,  provided  both  knife  edges 
lie  along  the  line  between  antennas  for  which  the  initial  calculation  would  otherwise 
yields  a  loss  associated  with  two  attenuation  factors  of  (1/2)  or 
loss  =  -20  logj  o  [(l  /  2)(l  /  2)]  =  \2dB . 
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RANGE  COORDINATE  (M) 

Figure  28. 


RANGE  COORDINATE  (M) 

Figure  29. 
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The  multiplicity  of  knife  edge  for  this  particular  simulated  profile  is  also  somewhat  due 
to  the  choice  to  model  a  structure  with  artificial  polygonal  segments  with  exact  linear 
height  functions  as  a  function  of  the  range.  Since  the  coordinates  understood  to  be 
elevations  above  sea  level  necessarily  exclude  earth  curvature,  the  construction  of  the 
actual  optical  obstructions  essentially  transforms  the  linearly  modeled  terrain  segment 
with  zero  curvature  into  convex  curves  with  a  negative  curvature  equivalent  to  that  of  the 
4/3  earth.  The  curvature  detection  algorithm  described  earlier  correctly  detects  the 
corresponding  portions  of  terrain  with  continuous  negative  curvature  and  the  associated 
terrain  points  are  processed  as  knife  edges.  Since  the  terrain  profiles  are  elevations  above 
sea  level,  optically  linear  constructions  of  zero  curvature  can  only  be  constructed  by 
including  a  compensating  positive  curvature  directly  into  the  terrain  profile.  If  the  terrain 
spacing  is  increased  by  a  factor  of  two  to  the  more  typical  spacing  approaching  15 
seconds  and  if  the  terrain  profile  is  adjusted  as  suggested  here  to  model  surfaces  of  zero 
curvature  in  the  optical  environment  of  a  4/3  earth  atmosphere,  the  plot  appearing  in 
Figure  29  results.  The  less  lossy  nature  of  the  resulting  propagation  modeling  is  readily 
apparent  from  the  spacing  of  contour  lines  which  is  2dB  in  both  Figure  28  and  Figure  29. 
In  view  of  its  effectiveness  in  the  reproduction  of  exact  results  for  two  knife  edges 
covering  a  wide  range  of  knife  edge  proportions,  it  is  realistic  to  expect  the  gain  factor  of 
Deygout  to  critically  assist  in  providing  results  similar  to  that  of  Figure  29  without  the 
uncertainties  of  artificial  terrain  or  knife  edge  spacing  restrictions.  For  truly  linear  terrain 
obstacles  as  intended  in  this  particular  test  profile,  modeling  by  perfectly  absorbing  knife 
edges  is  less  realistic  and  it  becomes  necessary  to  introduce  the  methods  of  GTD, 
discussed  in  the  following  section,  to  properly  model  the  associated  diffraction. 

The  juxtaposition  of  linearly  uniform  segments  defines  a  diffraction  center  characterized 
by  an  interior  wedge  angle,  a  geodesic  orientation  with  respect  to  an  effective  transmitter 
source  and  therefore  admits  to  a  modeling  consistent  with  the  wedge  diffraction 
formalism  of  Koyoumjian  and  Pathak  [6].  A  multiplicity  of  such  constructs  justifies  the 
standard  serialization  of  multiple  wedge  diffraction  centers  inherent  in  the  methods  GTD 
in  which  effective  serialized  source  geometries  are  defined  and  caustic  flux  preserving 
spreading  factors  are  well  established.  The  existence  of  terrain  realistic  impedance 
boundary  conditions,  serialized  geometries  which  introduce  overlapping  shadow 
boundaries  and  the  presence  of  horizontal  polarization  especially  in  association  with 
metallic  surfaces  requires  a  supplementing  of  GTD  with  slope  diffraction  and  its  semi¬ 
heuristic  adaptations  to  finite  ground  constants  [3],  [4],  [7],  The  desirability  of  GTD  was 
discussed  in  an  earlier  related  progress  report  including  coverage  plots  of  loss  in  the 
vertical  propagation  plane  for  a  truncated  wedge  Figure  17  and  for  the  rounded  wedge 
appearing  in  Figure  28  and  Figure  29.  The  functional  dependence  of  the  associated  GTD 
formalisms  on  source  coordinates  (ray  path  magnitude  and  angle  of  approach)  and 
corresponding  field  point  coordinates  eliminates  the  uncertainties  of  knife  edge  spacing 
tolerance  for  such  obstacles  exhibiting  the  requisite  polygonal  geometry.  Results  for  the 
same  calculation  at  300  MHz  (high  conductivity;  horizontal  polarization)  with  full  access 
to  the  methods  of  GTD  and  slope  diffraction  available  to  the  TMEM  GTD  Addendum 
appear  in  the  upper  left  panel  Figure  30  (and  for  5000  MHz  and  in  Figure  31)  The 
accompanying  panels  of  Figure  30  display  results  for  the  same  problem  but  using 
alternate  propagation  models,  TIREM,  APM  and  parabolic  equations  (VTRPE).  The 
TMEM  GTD  addendum  appears  to  provide  an  adequate  approximation  of  the 
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numerically  exact  results  shown  in  VTRPE.  Since  the  underlying  optical  geometrical 
parameters  for  computational  GTD  as  dictated  by  inter-relational  geodesic  geometry  are 
direct  generalizations  of  those  required  for  the  basic  diffraction  elements  presently 
implemented  in  the  developing  proposed  model,  the  extension  to  GTD  and  higher  order 
adaptations  is  expected  to  be  a  minimally  invasive  enhancement  of  current  methods,  fully 
implemented  within  the  established  computational  infrastructure. 

Finally  the  existence  of  inhomogeneous  earth  ground  constants  as  most  conspicuously 
exemplified  by  transitions  between  land  and  sea  will  require  an  additional  flexibility  to 
appropriately  adjust  ground  constants  and  weighted  partial  paths  loss  contributions  to  be 
in  conformity  with  established  limiting  quantitative  characteristics  of  path  loss  for 
propagation  over  mixed  land  and  sea  paths. 
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Figure  31.  5  GHz 

5.2.  Performance  Enhancements 

This  phase  of  research  and  development  has  focused  on  the  development  of  a  physically 
accurate  propagation  model.  As  the  model  is  slated  for  use  in  a  real-time  system,  it  must 
run  as  such.  Thus,  significant  improvements  to  the  speed  and  efficiency  of  the  algorithm 
are  planned. 
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6.  Schedule 


8.  Personnel  Support 

The  Principle  Investigator  for  this  project  is  Dr.  Stephen  A.  Fast.  Serving  as  Research 
Scientist  leading  the  Government  Service  Division  of  Remcom  Inc.  is  Stephen  A.  Fast. 
Dr.  Fast,  who  holds  a  Ph.D.  in  mathematics,  has  previously  occupied  Research  Associate 
positions  at  the  Applied  Research  Laboratory  at  The  Pennsylvania  State  University,  the 
Applied  Research  Laboratory  at  the  University  of  Texas  at  Austin,  the  Atlantic  Research 
Corporation’s  Electromagnetic  Environmental  Test  Facility,  and  the  U.S.  Army  Signal 
Corps.  The  work  done  at  Remcom  has  thoroughly  benefited  from  Dr.  Fast’s  academic 
background  and  expertise.  His  research  pursuits  include:  electromagnetic  propagation 
modeling  and  techniques,  cellular  phone  modeling  methods,  evaluation  of  tactical 
communications  and  intelligence  systems.  These  intellectual  enterprises  have  shaped  a 
remarkably  broad  range  of  professional  accomplishments  at  Remcom,  such  as  the 
development  of  the  Electromagnetic  Propagation  Integrated  Resource  Environment 
(EMPIRE),  and  a  single  point  urban  geolocation  method  based  on  urban  propagation 
predictions.  Dr.  Stephen  Fast  has  also  led  the  Remcom  in  a  highly  successful 
collaboration  with  the  PSU  Meteorology  Department,  which  focused  on  the  development 
of  atmospheric  interpolation  methods  for  anomalous  conditions  that  may  be  used  in 
support  of  radio  wave  propagation. 

Supporting  Dr.  Fast  is  Dr.  C.  Clay  Marston.  Dr.  Marston  recently  joined  Remcom  after 
serving  as  the  primary  theoretical,  algorithmic,  and  software  developer  of  TIREM  at 
Alion  Science  and  Technology  Corporation.  Dr.  Marston  holds  a  Ph.D.  in  physics  and 
has  conducted  fundamental  research  in  quantum  mechanics  and  theoretical 
electromagnetics.  His  current  research  interests  include  numerical  methods  in 
electromagnetics  especially  diffraction  and  scattering.  His  work  improving  the  methods 
of  TIREM  and  other  ray  based  propagation  models  have  resulted  in  new  techniques  for 
more  rigorous  methods  that  are  computationally  efficient. 

M.  Jeff  Barney  is  an  undergraduate  electrical  engineering  student  at  The  Pennsylvania 
State  University.  Jeff  is  interested  in  terrestrial  propagation  and  has  investigated  the 
various  differences  of  many  RF  propagation  models  including  VTRPE,  APM,  ITM, 
TIREM,  Prophet,  RPO,  and  MMwave.  Additionally,  Jeff  has  developed  methods  to 
derive  worldwide  foliage  maps  used  to  support  RF  propagation  through  foliage. 

9.  Publications 

The  results  of  the  first  portion  of  this  work  are  being  prepared  for  submission  to  Radio 
Science. 

10.  Interactions/Transitions 

Soon  after  the  end  of  the  first  year,  the  results  of  the  first  phase  shall  be  transitioned  to 
the  JCS. 
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11.  New  discoveries,  inventions,  or  patent  disclosures 

N/A 

12.  Honors/Awards 

N/A 
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13.  Appendix  -  Earlier  Progress  Report 

13.1.  Accomplishments/New  Findings 

13.1.1.  Baseline  run  times  and  model  comparisons 

Prior  to  implementing  any  new  methods  or  algorithms,  existing  models’  performance  was 
investigated.  Several  benchmarks  test  cases  were  run  to  compare  run  times  and 
propagation  predictions.  Thus,  those  areas  requiring  improvement  were  identified.  Thus 
far,  we  have  explored  the  terrain  extraction  performance,  compared  the  predictions  for 
various  models  for  a  few  standard  test  cases,  and  made  performance  comparisons  using 
the  most  promising  existing  models  and  algorithms  to  apply  to  the  JCS. 

13.1.1.1.  Terrain  extraction  study 

There  are  several  methods  to  compute  the  points  along  a  great  circle  path.  To  extract 
terrain  between  two  points,  the  range  and  bearing  between  them  must  be  computed.  Then 
along  the  great  circle,  latitudes  and  longitudes  must  be  found  and  corresponding  terrain 
computed  from  the  National  Imagery  and  Mapping  Agency  (NIMA)  Digital  Terrain 
Elevation  Data  (DTED).  Below  we  discuss  the  methods  that  could  be  applied.  In 
accordance  with  the  JCS  requirements,  only  1  km  DTED  terrain  data  will  be  used  and  the 
maximum  area  over  which  terrain  will  be  used  in  400x400  square  kilometers. 

13.1.1.1.1.  Modified  Rainsford  Method 

The  modified  Rainsford  method  is  one  of  the  most  accurate  methods.  This  method 
utilizes  a  real  earth  model,  i.e.  WGS-84,  and  an  iterative  algorithm  to  provide  very 
accurate  calculations  for  the  forward  and  inverse  problems,  i.e.  computing  a  point  on  the 
earth  given  the  range  and  bearing  from  a  given  point  and  computing  range  and  bearing 
given  two  points  on  the  earth.  However,  with  respect  to  performance  the  iterative  nature 
of  the  algorithm  can  cause  unnecessary  increases  in  time  and  complexity,  especially 
when  a  tolerable  level  of  error  can  be  achieved  with  less  complex  algorithms. 

13.1.1.1.2.  Linear  Interpolation 

This  method  simply  treats  the  latitude  longitude  grid  as  a  Cartesian  coordinate  system 
with  latitude  representing  X  and  longitude  representing  Y.  Using  this  method  the  end 
points  of  a  bearing  act  as  the  end  points  of  a  line  between  the  two  locations  and  using  the 
parametric  representation  of  this  line  we  can  find  the  location  of  any  points  in  between. 
This  method  is  simple  and  very  fast;  however  it  can  produce  large  errors  as  distance  from 
the  equator  increases. 

13.1.1.1.3.  Great  Circle 

This  method  uses  a  perfect  sphere  to  represent  the  surface  of  the  earth.  This  is  a  more 
reasonable  approximation  than  the  flat  earth  model  used  in  the  Linear  Interpolation 
method,  but  is  not  as  accurate  as  the  real  earth  model  used  in  the  modified  Rainsford 
method.  Using  this  method  is  a  compromise  between  the  two  previous  methods 
providing  more  accuracy  with  a  slight  performance  penalty.  This  method  primarily  uses 
trigonometric  functions. 
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13.1.1.1.4.  Test  Results 

To  determine  the  potential  run  time  improvement  over  the  Modified  Rainsford  Method  of 
the  other  methods,  a  test  was  run  in  Nellis  AFB  and  surrounding  area.  The  test  was 
designed  to  approximate  operation  in  a  400  km  X  400  km  area.  For  simplicity  we  chose 
a  slightly  larger  area  that  spanned  4  degrees  X  4  degrees.  This  area  was  sampled  on  the 
diagonal  for  a  total  of  100000  points.  The  results  are  depicted  in  Table  1  below. 


Table  1 


MRM  GC 

Linear 

Max  Difference  in  0  m  181.781m 

meters  from  Modified 
|  Rainsford  Method  j 

4957  m 

Speed  up  compared  to'  N/A  5/3 

MRM  ! 

o 

13.1.1.1.5.  Optimization 

Further  optimization  may  be  performed  by  allocating  fewer  temporary  variables  and 
precomputing  certain  values  for  the  Great  Circle  method.  Another  possible  solution  is  to 
combine  the  methods  in  a  dynamic  way  such  that  depending  on  the  area  of  the  world  and 
path  length.  A  combination  between  methods  may  allow  a  better  balance  of  performance 
and  accuracy.  Investigation  of  optimizations  for  the  Modified  Rainsford  Method  has  not 
been  completed  and  a  closer  investigation  of  the  method  may  reveal  some  possible 
optimizations  such  as  pre  calculation  of  sin/cos  tables  and  using  them  as  a  lookup  for  in 
subsequent  queries. 

13.1.1.1.6.  Recommendations  and  Conclusions 

Although  the  linear  interpolation  is  significantly  faster,  the  error  will  be  almost  5  terrain 
points  at  long  ranges.  The  great  circle  method  would  interpolate  between  the  same  terrain 
points  and  resulting  in  an  acceptable  error.  Thus  we  have  begun  building  our  software 
using  a  hybrid  method  between  linear  interpolation  the  great  circle  method. 

13.1.1.2.  Validation  comparisons 

Although  there  have  been  many  papers  comparing  various  propagation  models  and 
algorithms,  e.g.  knife-edge  and  other  diffraction  methods,  very  few  if  any  quantitative 
numerical  comparisons  have  been  done.  Thus  there  is  no  hard  data  on  which  to  select  a 
propagation  modeling  algorithm.  To  assist  selecting  a  proper  method  for  the  JCS,  a  few 
simple  test  cases  were  selected  and  some  of  the  currently  available  propagation  models 
were  compared  using  these  cases.  Before  proceeding  to  the  comparisons,  the  various 
methods  and  algorithms  for  each  propagation  model  are  discussed  below. 

Table  2  gives  a  brief  summary  of  the  models  and  their  capabilities. 
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Model  I  Frequency  |  Ground  |  Sky  wave  I  Atmosphere  I  Ground  I  Version  I  Algorithm  I  Comments 
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13.1.1.2.1.  TIREM 

The  Terrain  Integrated  Rough  Earth  Model  (TIREM)  is  a  general  purpose  point  to  point 
model  applicable  from  2  MHz  to  20  GHz  developed  by  the  Joint  Spectrum  Center  (JSC). 
TIREM  models  the  ground  wave,  terrain  reflection  and  diffraction,  and  troposcatter  in 
benign  atmospheres.  TIREM  computes  the  paths  between  the  terminal  points  of  a  radio 
link  and  applies  the  correct  theoretical  model  depending  on  the  situation.  When 
computing  rays,  the  effective  earth’s  radius  is  required  and  thus  TIREM  uses  the  surface 
refractivity  to  compute  the  effective  earth’s  radius.  Knowing  the  ray  paths  TIREM 
computes,  rain  and  foliage  attenuation  can  be  computed.  TIREM  is  one  of  the  most 
thoroughly  validated  models.  The  TIREM  validation  report  is  available  from  JSC,  the 
TIREM  sponsor,  which  outlines  over  10,000  real  world  measurements  used  to  validate 
TIREM.  According  to  the  report,  the  mean  error  in  TIREM  vs.  measured  data  is  less  than 
1  dB  with  standard  deviation  around  10  dB.  Further,  the  TIREM  developers  have  built 
TIREM  to  run  fast. 

TIREM  employs  a  hybrid  method:  switching  between  various  algorithms  as  the 
conditions  merit.  In  regions  where  more  than  one  propagation  mode  may  be  present, 
TIREM  selects  the  mode  with  the  least  loss. 

TIREM’s  algorithms  are  very  dependent  on  correct  terrain  sampling  resolution.  TIREM 
selects  the  radio  horizon  and  the  diffraction  and  reflection  point  by  searching  the  terrain 
points.  Thus,  JSC  recommends  that  TIREM  is  given  terrain  elevation  points  that  are  no 
more  the  250  meters  apart.  Unless  prominent  terrain  features  have  been  identified  and 
included  in  terrain  profiles,  JSC  recommends  3  arc  second  DTED  be  used.  However, 
terrain  reflection  is  only  considered  if  the  Fresnel  zone  intersects  the  terrain.  The  Fresnel 
zone  is  the  ellipse  of  points  such  that  the  path  difference  between  the  direct  path  and  any 
reflection  path  with  a  reflection  point  within  the  ellipse  is  less  than  Vi  a  wavelength.  For 
practical  purposes,  this  is  why  TIREM  results  do  not  show  the  strong  lobing  stucture.  The 
TIREM  Fresnel  zone  approximation  for  reflection  points  means  that  only  the  structure 
present  below  the  peak  of  the  first  lobe  is  modeled. 

For  terrain  diffraction,  TIREM  applies  the  Epstein-Peterson  knife-edge  diffraction 
method.  Depending  on  the  frequency  and  the  knife-edge  height  and  sharpness,  knife- 
edge  diffraction  may  underestimate  the  loss  in  the  shadow  region.  TIREM  compensates 
for  this  by  using  a  hybrid  algorithm  that  alleviates  this  problem  in  most  cases.  Thus,  the 
error  in  TIREM  has  been  reduced  to  only  a  few  dB  in  areas  in  the  shadow  well  above  the 
terrain. 

Often  in  deep  shadow  regions  the  troposcatter  path  will  have  less  loss.  TIREM  uses  the 
CCIR/ITU  method  to  predict  troposcatter  loss.  This  empirical  method  gives  only  mean 
results,  and,  thus,  troposcatter  results  may  vary  significantly  depending  on  the  area  of  the 
world  and  time  of  year.  In  general,  TIREM’s  method  will  produce  reliable  answers  that 
communication  planners  can  use  to  design  troposcatter  links. 

The  most  notable  TIREM  anomaly  is  found  in  the  shadow  region  of  a  broad  hill.  When 
computing  the  losses  for  a  broad  hill,  TIREM  often  identifies  more  than  one  knife-edge 
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on  the  hill.  This  causes  odd  shadowing  behavior  between  the  shadow  cast  by  the  first  and 
second  knife-edges. 

13.1.1.2.2.  ITM 

The  Institute  of  Telecommunication  Sciences  (ITS)  Irregular  Terrain  Model  (ITM)  is 
maintained  by  George  Hufford.  ITM  uses  ray  methods  together  with  parametric 
formulations  to  compute  the  loss  due  to  reflection,  diffraction,  and  troposcatter.  ITM  is 
applicable  between  20  MHz  and  20  GHz.  ITM  has  been  widely  accepted  in  the  radio  and 
television  broadcast  industry.  This  is  due  in  large  part  to  the  numerous  measurements 
taken  by  ITS  and  used  to  validate  the  model  and  develop  broadcast  statistical  parameters 
such  as  reliability  and  confidence  of  a  radio  link. 

The  reference  attenuation  is  not  a  smooth  function  of  distance,  but  is  divided  into  three 
regions.  The  line-of-sight  region  extends  to  the  smooth-earth  horizon  distance.  Here,  the 
bulge  of  the  earth  does  not  interrupt  the  direct  radio  waves  and  the  predicted  attenuation 
is  based  on  two-ray  reflection  theory  and  extrapolated  diffraction  theory.  In  the  next 
region,  the  predicted  attenuation  is  a  weighted  average  of  knife-edge  and  smooth  earth 
diffraction  calculations.  The  weighting  factor  is  a  function  of  frequency,  terrain 
irregularity,  and  antenna  heights.  Far  outside  the  radio  horizon,  forward  troposcatter 
effects  dominate.  An  adjustment  is  made  to  the  reference  attenuation  for  climate  to 
convert  the  reference  attenuation  to  an  all-year  median  attenuation. 

For  the  diffraction  region  only  the  transmitter  and  receiver’s  horizons  are  used.  The  loss 
is  computed  using  a  weighted  average  between  no  more  than  two  knife  edges  and 
Vogler’s  smooth  earth  methods.  The  weighting  factor  is  an  empirically  determined 
function  of  frequency,  the  antenna  heights,  and  terrain  parameters. 

Because  the  diffraction  rays  ignore  the  reflected  component,  the  deep  shadow  shows 
significantly  less  loss  than  is  actually  present.  Further,  as  only  the  endpoints’  horizons  are 
applied,  diffraction  effects  from  more  than  two  terrain  obstacles  will  not  be  evident  in 
ITM  predictions. 


13.1.1.2.3.  APM 

The  Advance  Propagation  Model  (APM)  was  developed  by  Amalia  Barrios  from 
SPAWAR  SSC-SD.  APM  uses  a  hybrid  parabolic  equation  (PE)  method  to  compute  the 
attenuation  from  100  MHz  to  20  GHz.  APM  models  terrain  reflection  and  diffraction, 
atmospheric  ducting,  and  troposcatter.  APM  uses  a  flat  earth  (FE)  ray  tracing  algorithm 
near  the  transmitter.  When  the  flat  earth  assumption  is  no  longer  valid,  APM  uses  a  ray 
optics  (RO)  method,  i.e.  ray  tracing  over  a  curved  earth.  When  the  transmitter  is  over 
rough  terrain,  neither  the  FE  nor  RO  methods  are  applied.  When  the  FE  and  RO  methods 
are  no  longer  valid,  APM  applies  the  parabolic  equation  (PE)  near  the  ground  and  an 
extended  optics  (XO)  method  above  the  PE  region. 

To  understand  the  importance  of  APM’s  hybrid  method,  first  consider  how  the  parabolic 
equation  method  works.  In  PE  models,  the  transmitter  antenna  pattern  is  converted  to  a 
field  distribution  along  a  vertical  line  located  at  the  calculation  origin,  i.e.  the  transmitter 
or  zero  range.  The  initial  field  is  then  marched  (stepped)  in  range  using  a  Fourier  method 
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and  an  atmospheric  correction.  Accordingly,  APM  and  other  PE  models  naturally  fill  a 
vertical  plane  with  loss  data  as  they  compute.  This  means  that  PE  models  are  not  efficient 
when  computing  an  isolated  point-to-point  calculation,  but  are  excellent  for  computing 
losses  for  a  large  numbers  of  height  points.  Within  each  range  step,  the  field  is  Fourier 
transformed  forward,  a  refractivity  correction  applied,  and  transformed  back  to  the 
original  domain.  For  the  Fourier  method  to  be  accurate,  it  must  have  a  height  grid 
resolution  of  around  0.1  X.  Since  the  vertical  extent  of  the  field  cannot  go  to  infinity,  a 
maximum  height  is  determined  in  the  algorithms,  and  methods  are  applied  to  eliminate 
effects  of  the  upper  boundary  of  the  calculation  region.  The  maximum  calculation  height 
is  related  to  the  vertical  beamwidth  that  the  algorithm  can  support;  more  correctly  the 
minimum  and  maximum  angles  in  which  energy  is  allowed  to  propagate.  Needless  to  say, 
full  PE  methods  with  wide-angle  capability  and  small  range  steps  are  computationally 
expensive.  Thus,  APM’s  hybrid  method  allows  it  to  reduce  the  run  time  penalties 
inherent  to  PE  methods.  To  get  further  run  time  improvements,  APM  also  sacrifices  a 
wide-angle  capability  and  some  accuracy  resulting  from  using  larger  range  steps. 

The  APM  propagation  model’s  narrow  angle  assumption  results  in  null  data  for  elevation 
angles  not  within  specified  limits.  For  elevation  angles  larger  than  APM’s  maximum, 
free  space  loss  can  be  substituted.  However,  because  of  possible  errors  caused  by  terrain 
reflection,  it  is  difficult  to  compensate  for  null  APM  data  for  angles  below  the  minimum 
elevation  angle.  This  will  be  most  noticeable  when  the  transmitter  is  high  and  over  rough 
terrain. 

Using  APM,  shadows  behind  obstacles  often  have  losses  similar  to  those  observed  with 
knife-edge  diffraction,  possibly  because  of  the  angle  limitation  and  step  size  selection. 
Thus,  APM  underestimates  the  loss  in  the  deep  shadow.  APM  uses  a  fixed  step  size,  so, 
occasionally;  the  internal  APM  step  will  not  coincide  with  the  diffraction  edge.  Thus  the 
shadow  boundary  in  APM  is  different  than  that  in  models  which  ensure  the  diffraction 
point  is  identified  and  used. 

As  with  all  PEs,  the  initial  field  in  APM  is  found  from  the  transmitter’s  antenna  pattern. 
Thus,  unlike  TIREM,  the  transmitter’s  antenna  pattern  cannot  be  applied  after  the 
calculation.  Also,  because  of  the  nature  of  the  PE  algorithm,  APM  does  not  provide  data 
to  an  external  interface  in  order  to  allow  the  elevation  angle  dependence  of  the  receive 
antenna  to  be  considered.  Thus,  it  is  difficult  to  receive  antenna  patterns  when  using 
APM. 

The  APM  API  requires  that  the  calling  application  provide  the  refractive  data.  In  order  to 
interpolate  the  refractive  data  at  the  PE  range  steps,  APM  assumes  that  the  refractive  data 
can  be  interpolated  using  corresponding  indices  in  the  profile  arrays.  Therefore,  when 
using  APM,  the  calling  program  is  responsible  for  arranging  the  refractive  profiles 
accordingly,  so  that  the  interpolation  is  physically  correct.  In  addition  to  refractive  profile 
data,  APM  uses  the  humidity  and  temperature  to  compute  absorption. 

APM  uses  the  wind  speed  to  determine  the  surface  roughness  of  waves  over  water.  The 
wind  speed  is  related  to  the  wave  height  and  spacing,  and  consequently,  the  surface 
roughness.  The  surface  roughness  is  then  mapped  to  an  angle  dependent  reflection 
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coefficient.  Thus,  APM  essentially  accounts  for  surface  roughness  by  changing  the 
reflection  coefficient.  Thus,  energy  is  not  scattered,  but  the  effects  of  scatter  on  the 
forward  propagated  energy  is  considered.  APM  allows  the  wind  speed  to  be  set  as  a 
function  of  range.  Thus,  the  surface  roughness  can  be  varied  along  a  path.  As  well  as 
computing  the  surface  roughness  effects  over  water,  APM  is  equipped  to  take  advantage 
of  range  dependent  dielectric  constants. 

Lastly,  one  manifestation  of  APM’s  hybrid  scheme  is  an  occasional  inconsistency  in  the 
results  as  the  boundary  changes  from  the  PE  to  the  XO  regions.  In  unpredictable 
situations,  a  feature  or  structure  will  slightly  and  abruptly  change  direction  when  crossing 
the  PE/XO  boundary. 

13.1.1.2.4.  VTRPE 

The  Variable  Terrain  Radio  Parabolic  Equation  (VTRPE)  model  is  one  of  the  most 
versatile  and  precise  models  available.  VTRPE  was  developed  by  Frank  Ryan  from 
SPA  WAR  SSC-SD.  VTRPE  can  model  the  ionospheric  sky  wave,  terrain  reflection  and 
diffraction,  ducting,  and  the  ground  wave. 

VTRPE  is  a  full  PE,  i.e.  no  approximations  are  made  to  the  PE  method  and  no  hybrid 
methods  are  used.  VTRPE  uses  a  wide-angle  propagator,  i.e.  energy  at  very  large 
elevation  angles  is  propagated  (as  compared  with  APM  which  has  an  upper  elevation 
angle  limit).  VTRPE  examines  the  field  at  each  range  step  to  determine  the  Fourier 
transform  size  required  to  correctly  accomplish  the  computations.  Further,  VTRPE 
evaluates  the  fields  at  each  range  step  to  determine  the  maximum  acceptable  step  size.  By 
allowing  the  Fourier  transform  size  and  the  range  step  to  be  variable,  VTRPE  keeps  the 
results  accurate  and  optimizes  the  model  execution  time.  Thus,  if  a  transmitting  antenna 
has  high  elevation  angle  contributions,  e.g.  an  isotropic  antenna,  the  Fourier  transform 
size  will  be  large  and  the  range  step  size  small  causing  longer  than  potentially  required 
run  times. 

As  with  APM,  VTRPE  can  model  the  effects  of  range  dependent  atmospheres  and 
dielectric  constants.  VTRPE  can  either  call  out  to  an  external  program  for  the  refractive 
index  data  at  each  range  step  or  have  multiple  profiles  passed  to  it.  In  comparison  to  the 
atmosphere,  VTRPE  models  the  ionosphere  internally.  The  model  VTRPE  uses  is  the 
FAIM-D  ionosphere  model  for  the  E  and  F  layers  and  the  smoothly  transitions  to  the  IRI 
model  to  include  the  D  Layer. 

Antenna  patterns  for  VTRPE  are  treated  similarly  to  APM,  that  is,  the  transmit  antenna 
pattern  is  applied  to  compute  the  initial  field.  However,  since  VTRPE  can  model  the 
ground  wave,  VTRPE  models  the  complex  antenna  ground  interaction  required  to 
properly  excite  the  ground  wave  for  arbitrary  antennas.  Thus  VTRPE  requires  that  for 
antennas  model  the  ground  plane  is  not  included  in  the  pattern.  Additionally,  for  non  HF 
receive  antennas,  since  VTRPE  only  offers  loss  and  no  angular  distribution  of  the  energy 
at  the  receive  point. 
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13.1.1.2.5.  Comparisons 

To  understand  the  differences  between  the  models,  some  simple  test  cases  using  terrain 
that  allows  the  differences  to  be  isolated  were  used.  The  configurations  for  these  tests  are 
described  below.  In  all  cases  the  frequency  chosen  for  the  comparisons  was  300  MHz 
with  vertical  polarization.  APM,  VTRPE,  TIREM  and  ITM  were  compared. 

Two  wedges  were  chosen  to  better  quantify  the  effects  of  a  second  diffraction.  The  first 
double  wedge  case  uses  two  wedges  with  difference  maximum  elevations.  This  allows 
the  testers  to  see  if  the  shadow  boundary  caused  by  diffraction  from  the  second  wedge  is 
properly  modeled.  The  transmitter  height  is  10  meters,  the  min  range  =  5  km,  max  range 
=  50  kin,  min  height  =  1  meter,  max  height  =  1,001  meters  MSL,  number  of  ranges  = 
181,  number  of  heights  =101,  and  the  terrain  profile  in  range  in  km  and  height  in  meters 
is  (0,  0),  (6.25,  200),  (12.5, 0),  (25, 0),  (31.25,  150),  (37.5,  0),  and  (50,  0). 

The  test  results  for  the  two  uneven  wedges  are  show  in  Figure  25.  First  notice  the  VTRPE 
results.  The  interference  pattern  from  the  direct  and  reflected  rays  is  clearly  present  in  the 
illuminated  region.  Also,  we  see  the  diffracted  and  reflected  ray  interference  in  the 
shadow  regions.  Evidence  of  reflected  energy  off  of  the  front  face  of  the  second  wedge  is 
also  seen.  APM  does  almost  as  nice  a  job  as  VTRPE.  With  APM  the  interference  pattern 
in  the  shadow  of  the  first  wedge  is  noticeably  bent  in  a  non-physically  justifiable  way. 
Also,  since  APM  uses  a  narrow  band  PE  propagator,  the  strong  lobing  found  before  the 
first  hill  is  not  present.  In  short  the  parabolic  equation  models’  results  demonstrate  why 
the  PE  method  is  widely  accepted  as  the  standard.  Only  high  fidelity  ray  methods  such  as 
those  used  in  GTD  Estimation  of  Loss  due  to  Terrain  Interaction  (GELTI)  will  give 
similar  results  as  the  PE  methods.  Thus  we  should  not  expect  TIREM  or  ITM  to  show  the 
same  features  as  VTRPE  and  APM.  First,  neither  TIREM  nor  ITM  model  the  reflected 
ray  (except  very  near  the  ground)  and  thus  the  lobing  structure  is  not  found  in  these 
results.  The  transition  to  the  shadow  boundary  is  significantly  sharper  in  these  models.  In 
TIREM  diffraction  is  applied  just  above  the  shadow  boundary  so  the  transition  is  more 
realistically  gradual  in  its  results.  ITM,  however,  does  not  model  this  at  all  and  the 
transition  is  very  abrupt.  Lastly,  in  the  shadow  of  the  first  hill,  TIREM  underestimates 
the  loss. 

The  second  double  wedge  case  (Figure  26)  uses  two  wedges  with  the  same  maximum 
elevations.  Again,  the  transmitter  is  10  meters  above  the  ground  and  min  range  =  5  km, 
max  range  =  50  km,  min  height  =  1  meter,  max  height  =  1,001  meters  MSL,  number  of 
ranges  =181,  number  of  heights  =101,  and  the  terrain  profile  in  range  in  km  and  height 
in  meters  is  (0,  0),  (6.25, 0),  (12.5, 450),  (18.75, 0),  (31.25, 0),  (37.5, 450),  (43.75, 0),  and 
(50,  0). 

For  the  two  even  wedges  (now  slightly  higher  and  further  from  the  transmitter)  we  notice 
some  interesting  differences  in  the  models’  behavior.  Again,  VTRPE  gives  results  that 
demonstrate  the  presence  of  all  of  the  correct  physics  are  found  in  the  algorithm. 
Interestingly,  in  the  VTRPE  results  we  see  the  energy  reflected  from  the  front  of  the  first 
wedge  and  the  diffraction  pattern  above  the  second  wedge;  features  not  seen  in  the  other 
models.  Unlike  the  first  case,  TIREM  shows  the  deep  shadow  expected  behind  the  first 
wedge.  However,  TIREM  shows  an  anomaly  emanating  from  the  second  wedge  near  the 
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shadow  boundary.  ITM  also  shows  an  unusual  feature  near  the  first  shadow  boundary. 
Further,  the  deep  shadow  behind  the  first  wedge  is  not  comparable  with  the  PE  results. 

Figures  27  and  28  show  the  results  of  the  two  remaining  test  cases.  The  first  shows  the 
interaction  with  a  rounded  hill.  As  before,  the  transmitter  is  10  meters  above  the  ground, 
the  min  range  =  5  km,  max  range  =  50  km,  min  height  =  1  meter,  max  height  =  1,001 
meters  MSL,  number  of  ranges  =  181,  number  of  heights  =  101,  and  the  terrain  profile  in 
range  in  km  and  height  in  meters  is  (0,  0),  (18.75,  0),  (20.312,  210),  (21.875,  320), 
(23.4375,  375),  (25,  390),  (26.5625,  375),  (28.125,  320),  (31.25,  90),  (32.8125,  0),  and 
(50,  0). 

This  test  case  challenges  the  models  to  represent  a  large  rounded  hill.  Often  this  test  will 
cause  a  model’s  diffraction  search  algorithms  to  give  erroneous  results.  In  fact,  this  test 
case  demonstrates  that  TIREM  may  be  unable  to  properly  select  the  diffraction  edges. 
Also  because  of  the  rounded  transition  to  the  shadow  beyond  the  wedge,  models  that  treat 
diffraction  as  a  pure  knife  edge  fail  to  give  dark  enough  shadows  near  the  ground  behind 
the  hill.  Thus  this  test  case  illuminates  ITM’s  method  for  computing  diffraction  and  the 
resulting  30-40  dB  difference  between  ITM  and  the  PE  models  in  the  deep  shadow.  As 
has  been  observed  in  the  previous  test  results,  the  PE  models  give  physically  reasonable 
results  for  this  test  case.  Further  VTRPE  again  models  all  of  the  propagation  features 
expected  by  the  physics  governing  this  situation. 

The  final  test  case  seeks  to  isolate  diffraction  from  the  other  phenomena.  In  this  case  a 
transmitter  at  1,278  meters  MSL,  vertical  polarization,  min  range  =  25  km,  max  range  = 
50  km,  min  height  =  1,248  meters,  max  height  =  3,248  meters  for  AGL  and  MSL, 
number  of  ranges  =  201,  number  of  heights  =  101,  the  terrain  profile  in  range  in  km  and 
height  in  meters  to  (0, 0),  (25,  1,  748),  and  (50, 0). 

After  reviewing  the  results  from  the  previous  test  cases,  the  results  for  the  final  test  are 
expected.  Again  the  PE  methods  demonstrate  the  underlying  physics  very  well.  TIREM 
matches  the  PE  methods  quite  well  in  the  shadow  region  while  ITM  underestimates  the 
loss  in  the  deep  shadow. 

These  test  results  confirm  that  correct  modeling  of  the  shadow  boundary,  reflection  lobes, 
and  deep  shadows  is  required  for  JCS  scenarios.  After  reviewing  these  test  results  and 
considering  them  together  with  the  platform  speed  and  environmental  parameter 
resolution,  it  is  clear  that  modeling  the  following  features  is  less  important  than  modeling 
the  aforementioned  first  order  features:  small  scale  lobing  in  the  illuminated  diffraction 
region,  smooth  accurate  transition  from  the  illuminated  to  the  shadow  region,  reflection 
structure  in  the  shadow  region,  and  reflection  from  the  front  faces  of  hills.  These  minor 
features  do  not  have  large  loss  variations  and  largely  occur  over  small  areas.  Thus  the 
time  resolution  of  the  simulation  and  corresponding  spatial  resolution  render  the  second 
order  effects  inconsequential  to  JCS  simulations. 
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13.1.1.3.  Run  time  comparisons 

Using  terrain  data  from  Edwards  Air  Force  Base,  the  run  times  for  36000  TIREM  and 
ITM  computations  were  compared.  All  of  the  following  tests  were  run  on  an  Intel 
Pentium  4  (3.06  Ghz)  with  512MB  RAM  and  a  512  KB  Cache  using  the  Redhat  Linux 
8.0  operating  system.  From  the  Sec®  CPU  2000  benchmark  data,  chips  found  in  the  JCS 
have  about  40%  of  the  performance  of  the  Pentium  4  used  for  this  test. 
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Although  these  tests  were  not  completely  comprehensive  and  the  results  preliminary, 
clearly  ITM  is  much  faster  than  TIREM.  So  once  the  core  methods  of  these  algorithms 
can  be  extracted  and  the  blatant  intolerable  errors  corrected,  they  can  be  written  in  the  C 
programming  language.  The  methods  currently  under  development  addressing  this  are 
discussed  below. 

13.1.2.  Proposed  methods  and  algorithm 

13.1.2.1.  Derivation  and  Discussion 

This  section  will  outline  the  theoretical  motivation  and  practical  justification  for  a 
number  of  features  presently  installed  in  the  developing  model.  The  first  subsection  will 
describe  the  method  for  enforcing  continuity  across  the  shadow  boundaries  associated 
with  initial  obstruction  of  inter-visibility  and  with  higher  order  shadows  involving 
multiple  diffraction  centers.  Subsection  [13.1.2.1.2]  will  present  a  qualitative  description 
of  a  method  employed  to  model  specular  reflections  accompanying  smooth  earth 
attenuation  of  initially  knife-edge  diffracted  fields.  Subsection  [13.1.2.1.3]  will  examine 
the  multiplicity  of  theoretical  and  practical  aspects  of  GTD  and  its  accepted 
generalizations  which  contribute  to  its  particular  applicability  to  the  present  model.  The 
final  section  will  detail  two  representative  auxiliary  algorithmic  constructs  which  have 
proven  vital  to  the  inevitable  realization  of  the  current  objectives. 

13.1.2.1.1.  Initial  and  Multiple  Knife  Edge  Shadow  Boundary  Continuity 

A  desirable  feature  of  the  proposed  model  is  a  continuous  interfacing  of  fields  initially 
propagated  over  a  smooth  earth  which  then  encounter  shadow  boundaries  introduced  by 
irregular  terrain  features  best  modeled  by  GTD  or  single  knife  edge  diffraction. 
Continuity  across  the  shadow  boundary  of  a  smooth  earth  diffracted  field  attenuated  by  a 
knife  edge  can  be  enforced  by  a  method  similar  to  that  found  effective  in  enforcing 
continuity  across  the  (second)  shadow  boundary  of  a  double  knife  edge  diffraction  as 
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modeled  by  Epstein  and  Peterson  [1].  Epstein  and  Peterson  [1]  model  the  double 
diffraction  regions  below  the  second  shadow  boundary  as  the  product  of  two  attenuations 
each  calculated  as  an  effective  single  knife  edge  attenuation  in  which  the  first  knife  edge 
serves  as  the  source  for  the  second.  Thus  in  the  Figure  29  curves  Rl23  and  /?223are 
associated  with  the  conventional  attenuation  factors  in  which  the  first  knife  edge 
contributes  the  constant  values  a0]2  along  curve  /?l23,  below  point  C,  and  the  second 

knife  edge  contributes  the  variable  attenuation  a12(v2)  along  curve  /?223  below  transition 
region  D. 

Above  the  second  shadow  boundary  the  total  attenuation  may  also  be  modeled  as  the 
product  of  two  terms,  the  first  of  which,  due  to  the  first  knife  edge  diffraction,  Ax  is  taken 

as  a01(v,)  along  curve  Rm  above  point  C,  with  a  lower  limit  set  by  its  values  at  the 
second  knife  edge  A,  =  a0l2 .  The  second  factor  due  to  the  diffraction  by  the  second  knife 
edge  A2  may  be  defined  as  one  for  values  of  the  clearance  ratio  exceeding  transition 
argument  v  =  0.78 ,  above  transition  region  D,  and  as  the  conventional  second  knife  edge 
attenuation  A2  =  an  (v2 )  for  arguments  less  than  0.78,  below  transition  region  D.  The 
activation  of  the  second  attenuation  at  .78  is  continuous  because  the  analytic  attenuation 
a12(v2  =.78)  is  identical  to  the  defined  value  of  unity  for  values  of  v2  exceeding  .78.  The 
loss  increases  continuously  as  a  function  of  v2  from  an  initial  loss  of  0  at  the  transition 
argument  of  v2  =  .78  to  values  of  loss  exceeding  0  for  values  of  v2  less  than  .78. 

Proceeding  similarly  for  the  shadow  boundary  associated  with  a  smooth  earth  diffracted 
field  attenuated  by  a  knife  edge,  the  first  attenuation  factor  is  assigned  its  conventional 
value  of  as0e  above  point  A  along  curve  Rm  and  is  assigned  the  limiting  (constant)  value 
of  <  below  point  A.  The  second  attenuation  factor,  due  to  the  first  knife  edge,  is  defined 
as  1  along  curve  Rm  for  arguments  exceeding  .78  (above  transition  region  B)  and  is 
otherwise  assigned  its  conventional  single  knife  edge  attenuation  a0l  along  curve  /?,  12 , 
below  transition  region  B.  The  success  of  this  mechanism  of  enforcement  of  continuity 
across  the  smooth  earth-knife  edge  shadow  boundaries  is  apparent  in  the  loss  contour 
plots  presented  in  this  report. 

13.1.2.1.2.  Derivation  of  Terrestrial  Specular  Knife  Edge  Reflection 
This  section  will  present  a  qualitative  description  of  a  method  employed  to  model 
specular  reflections  accompanying  smooth  earth  attenuation  of  initially  knife-edge 
diffracted  fields.  The  smooth  earth  field  propagation  formulation  of  Norton  [10]  involves 
the  summation  of  a  direct  ray  and  a  reflected  ray  which  is  corrected  for  finite  ground 
constants  by  a  surface  wave 
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Factoring  out  the  direct  ray,  this  may  be  written  explicitly  as  the  product  of  the  direct  ray 
and  a  smooth  earth  interference  attenuation  factor: 
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If  the  complex  direct  and  reflected  prefactors  in  Eq.(l)  are  replaced  by  smooth  earth- 
knife  edge  attenuated  field  emanating  from  the  first  knife  edge,  one  obtains 

E(se’ke)  (v)  =  Eke  (v, )  +  Eke  (v2  \DR  +  (1  -  R)f(P,  B)]  (3) 


where  quantities  within  the  bracketed  term  are  computed  relative  to  the  effective 
transmitter  provided  by  the  intervening  knife  edge,  and  where  the  Fresnel  integral 
representation  of  knife  edge  attenuation  assumes  the  form  [8] 
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The  Fresnel  integrals  appearing  in  this  more  general  problem  are  not  related  by  a 
multiplicative  phase  factor,  but  are  rather  explicit  functions  of  the  coordinates  associated 
with  the  real  and  virtual  receivers  introduced  by  image  theory  and  must  therefore  be 
evaluated  independently  as  in  Eq.(3).  The  above  approximation  is  implemented  by 
evaluating  the  knife  edge  attenuation  at  both  real  and  virtual  receivers,  evaluating  the 
bracketed  expression  in  Eq.(3)  with  arguments  consistent  with  a  transmitter  located  at  the 
intervening  knife  edge,  and  then  adding  the  resulting  generalized  direct  and  reflected 
components.  The  inclusion  of  the  surface  wave  component,  appearing  as  the  second  term 
in  Eq.  (3),  is  expected  to  eliminate  the  inordinately  high  losses  that  often  accompany  the 
mutual  cancellation  of  direct  and  reflected  fields  as  modeled  by  a  simple  two  ray  specular 
reflection  construction. 

13.1.2.1.3.  Justification  for  Selection  of  GTD  and  Extensions 

This  section  will  examine  the  appealing  features  inherent  in  GTD,  its  extension  to 
multiple  diffraction  and  its  reinterpretation  as  the  first  term  of  an  exact  expansion. 
Techniques  of  geometrical  optics  applicable  to  classical  diffraction  originated  by  Keller 
[5]  and  refined  by  Kouyoumjian  and  Pathak  [6]  and  practically  adapted  by  Luebbers  [7] 
and  Holm  [3] [4]  are  appealing  to  the  present  objective  and  foundational  to  the  example 
calculations  presented  in  this  report.  The  especially  appealing  features  of  GTD/UTD 
include  its  reduction  to  established  exact  first  order  approximations  for  both  single  and 
double  knife  edge  configurations  and  its  provision  of  a  multiple  diffraction  algorithm  for 
both  knife  edges  and  wedges  which  are  identical  to  established  first  order  terms  for  the 
more  general  multiple  diffraction  configurations. 

Kouyoumjian  and  Pathak  [6]  show  that  their  formulation  of  wedge  diffraction,  if  applied 
to  the  special  case  of  vanishing  interior  wedge  angle,  naturally  reduces  to  Sommerfeld's 
solution  for  diffraction  by  a  half  plane.  Also  in  the  context  of  applications  for  single  knife 
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edge  diffraction  and  as  applied  to  finite  ground  constants,  Holm  [4]  has  shown  that  a  self- 
consistent  modeling  of  the  four  UTD  partial  diffraction  coefficients  identifies  a  single 
term  as  being  associated  with  a  perfectly  absorbing  knife  edge.  This  result  is  also 
consistent  with  the  more  general  results  of  Sommerfeld  for  perfect  conductivity  of  either 
polarization. 

Extending  this  comparative  study  to  established  analytic  results  for  two  knife  edges,  in 
examining  Holm's  novel  expression  for  multiple  wedge  diffraction  as  a  sum  of  products 
of  one  dimensional  Fresnel  integrals,  it  becomes  apparent  that  for  the  special  case  of  two 
diffraction  centers  with  zero  interior  wedge  angle,  Holm's  expression  reduces  to 
Millington's  exact  derivation  [9]  of  the  correction  of  the  Epstein  Peterson  model  [1]  for 
attenuation  by  two  knife  edges.  In  particular  Millington  shows  that  the  Epstein  Peterson 

expression  is  exact  to  within  a  factor  of  sin  (y)  such  that  the  magnitude  of  the  doubly 
diffracted  field  is  written 


=  «012«12(V2)sinW  (5) 

where  a0l2  and  al2(v2)  are  attenuation  factors  associated  with  the  first  and  second  knife 

edges.  Millington's  essential  correction  of  the  Epstein  Peterson  formulation  is  the  pre¬ 
factor  multiplying  the  product  of  the  single  knife  edge  attenuations  which  in  his  notation 
is  written 


sin(y)= 


b{a  +  b  +  c) 

( a  +  b\b  +  c) 


ll/2 


(6) 


Converting  to  Holm's  notation  now  by  making  the  substitutions  a=  n,  b=r2  and  c=r3,  one 
obtains 


rM  +  r2  +  r3) 
(ri  +  r2\r2  +  r3 ) 


1/2 


(7) 


The  analogous  expression  derived  by  Holm,  in  Eq.(9)  of  Reference  [3],  intended  to  also 
express  double  knife  attenuation  as  a  product  of  single  knife  edge  attenuation  factors,  is 


Introducing  the  transformation  expression 


rir2 

n  +  r2 


^2  ~ 


r2+r3 


The  Holm  double  knife  attenuation  prefactor  becomes 

A  _\{n  +  r2  +  r3Xr,r2)(r2r3) 
+  r2)(r2  +  r3) 


which  reduces  to 


An 


(r,  +  r2  +  r3)r2 

(>i  +  r2)(r2  +  r,) 


1/2 


(9) 


(10) 


which  is  identical  to  Millington's  expression  ,  Eq.  (6). 

In  view  of  the  multiplicity  of  algorithms  [2]  proposed  for  modeling  multiple  knife  edge 
diffraction,  it  is  appealing  that  the  established  extension  of  GTD  to  multiple  diffraction 
has  been  identified  as  the  leading  term  of  Vogler's  exact  series  solution.  The  mechanical 
generalization  of  GTD  to  multiple  diffraction  when  reduced  to  wedges  with  zero  interior 
wedge  angle  provides  a  multiple  knife  edge  algorithm  with  all  the  inherent  advantages  of 
GTD.  Holm  [3]  has  shown  that  GTD  may  be  interpreted  as  the  first  term  resulting  from  a 
rearrangement  of  Vogler's  exact  multiple  knife  edge  expression  and  that  the  exact  result 
may  be  obtained  by  including  sufficiently  many  terms.  Holm  has  also  shown  that  the 
second  term  of  his  summation  corresponds  to  slope  diffraction  which  has  been 
established  as  the  leading  second  order  correction  of  first  order  GTD.  The  provision  by 
GTD  of  a  multiple  knife  edge  algorithm  recognizable  as  the  leading  term  of  established 
exact  expressions  for  multiple  diffractions  and  recognition  that  the  higher  order  terms 
themselves  consist  of  one-dimensional  integrals  reinforces  the  perception  of  its 
desirability  for  the  present  purposes.  Similar  remarks  apply  to  finite  wedge  angles  and 
therefore  an  appealing  continuity  of  formalism  is  realized  in  the  generalization  of 
multiple  knife  edge  to  multiple  wedge  diffraction.  Thus  GTD  serves  as  first  term  of  exact 
solution  expressed  as  a  summation  which  is  especially  relevant  to  previously  intractable 
problems  associated  with  overlapping  shadow  boundaries  in  that  they  may  be  solved  by 
the  simple  inclusion  of  more  terms. 

13.1.2.1.4.  Auxiliary  Constructions  Supportive  of  the  Current  Effort 

This  subsection  will  describe  examples  of  supporting  algorithms  developed  during  the 
time  frame  of  the  current  review  period.  The  first  portion  of  this  subsection  describes  the 
complexity  of  proper  computation  of  height  gain  functions  inherent  in  smooth  earth 
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propagation  which  has  required  numerical  realization  of  the  essential  methods  outline  by 
Norton  [10].  Finally,  a  method  for  eliminating  troublesome  interference  variability  is 
examined  at  the  conclusion  of  this  subsection. 

Smooth  earth  algorithm  of  Norton  with  complete  treatment  of  surface  wave,  height  gain 
functions,  and  derivative  conserving  intermediate  regime  interpolation  algorithm. 

To  assure  proper  treatment  of  those  aspects  of  terrestrial  propagation  best  modeled  as 
diffraction  over  a  smooth  earth,  the  formulations  of  Norton  [10]  have  been  adapted  in  the 
time  frame  covered  by  this  progress  report  to  computational  accessibility.  Norton's  form 
for  the  geometrical  optics  electric  field  propagated  over  a  smooth  planar  earth  is  given  by 

£=£H1+e“tts+<1"*>/(p’4  (11> 

where  A  is  the  phase  difference  between  direct  and  reflected  ray  paths  A  =  2 7i{r2  -  r,) /  / 1 

For  extended  planar  configurations  with  antennas  close  to  the  ground,  the  reflection 
coefficient  approaches  -1,  and  the  above  planar  earth  formulation  as  modeled  by  Norton 
approaches  an  intermediate  asymptotic  form  associated  with  complete  mutual 
interference  of  direct  and  reflected  waves  to  define  a  surface  wave  which  must  be  refined 
by  height  gain  factors: 


-ifcr, 

E(r)  -  y  2 f(p,b)z{^)\  02) 

where  z(qi);i  =  1 ,2  is  the  height  gain  factor  associated  with  the  transmitter/receiver  and 
where  qx  and  q2  are  scaled  antenna  heights  normalized  by  ground  constant  parameters. 
For  sufficiently  elevated  antennas,  the  argument  dependence  of  the  height  gain  factors 
must  be  generalized  to  include  another  argument,  the  ground  constant  and  frequency 
dependent  parameter  K  such  that  the  above  equation  is  written 

-/Ar(  _ 

£W= - — (13) 

For  sufficiently  elevated  antenna,  the  two  physical  form  of  the  diffracted  field  are  unified 
in  a  single  theoretical  formulation 

E(r)  =  ^-fl  +  e“  (I-  /(O/f/’.fi)]  (14) 

Ant\  r2l  JJ 
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where  primes  indicate  dependence  upon  spherical  earth  geometry  and  where  D  is  a 
divergence  factor  not  explicitly  present  in  the  planar  earth  formulation.  The  height  gain 
factors  in  the  extended  planar  surface  wave  and  in  the  more  general  spherical  earth 
formulation  may  be  identified  as  factors  that  appear  in  Bremmer's  exact  solution  which 
admits  to  a  Hankel  approximation  involving  solution  of  the  root  generating  equation 

tf,1/3(<f)  =  (3  if  ei(-bl2^ni)  H2x'\i)  (15) 

where  is  the  complex  Hankel  function  of  the  first  kind  of  real  order  v  and 

£=  J{-  2r0)  .  The  arguments  q,  and  K  of  the  height  gain  factors  z(qi ,  K) are  functions 

of  the  complex  quantity,  r0  appearing  in  the  above  equation  which  must  be  solved  to 
properly  compute  the  vital  height  gain  factors  in  Norton's  formulation.  The  above 
formulations  are  valid  within  a  restricted  range  of  the  horizon  dictated  by  the  limits  of 
geometrical  optics  such  that  the  proper  form  beyond  these  limits  is  computed  as  an 
interpolating  function  between  the  geometrical  optical  limits  and  the  final  diffracted  form 
approximated  by  the  first  terms  of  the  Bremmer  series.  Within  the  time  frame  covered  by 
this  progress  report,  the  graphical  methods  of  Norton  have  been  mathematically 
transformed  for  computational  accessibility,  the  proper  special  functions  for  computing 
the  surface  wave  have  been  collected  and  numerical  iterative  solutions  to  the  complex 
Hankel  root-generating  equation  have  been  developed  to  compute  the  complex  quantities 
upon  which  the  arguments  of  the  height  gain  functions  depend.  The  success  of  the 
implementations  is  apparent  from  the  replication  of  Figures  13-15  in  Norton's  paper. 

Elimination  of  numerical  precision  variability  associated  with  superposition  of  UTD 
direct  and  diffracted  components  with  essential  first  order  small  argument  phase 
difference  expression  at  overlapping  shadow  boundaries. 

Above  the  shadow  boundaries,  the  formalism  of  GTD  requires  that  the  direct  ray  serving 
as  source  to  the  diffraction  center  is  added  to  the  diffracted  ray  itself.  In  the  approach  to 
shadow  boundary  from  the  illuminated  half-plane,  the  diffracted  ray  approaches  a  form 
which  may  be  identified  as  -1/2  the  direct  ray  and  the  superposition  of  the  two  leads  to  a 
total  field  which  is(l/2)  times  the  direct  ray  for  the  characteristic  6  dB  loss.  Below  the 
shadow  boundary,  the  direct  ray  no  longer  contributes  and  the  diffraction  formalism 
imposes  a  change  in  sign  on  the  diffracted  ray  such  that  the  total  ray  is  purely  diffracted 
and  of  exactly  the  same  magnitude  as  the  total  field  infinitesimally  above  the  shadow 
boundary.  Coexisting  with  the  delicacies  of  field  continuity  across  the  shadow  boundary 
is  a  possibility  for  an  undesirable  numerical  variability  resulting  from  the  addition  of  the 
direct  and  diffracted  rays  which  is  attributable  to  small  phase  differences  of  the  two  wave 
forms.  The  interference  factor  of  the  direct  and  diffracted  components  can  be  represented 
by 


-ikrx 


Anr. 


[l  +(r!  p')DAe-lk[p'p'-r)] 


(16) 
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The  phase  factor  r~ik^p+p'^  essentially  involves  the  subtraction  of  one  quantity,  r,  the 
magnitude  of  the  direct  ray  from  the  sum  of  two  quantities  (/?+/?')  to  produce  a  phase 

difference  factor  which  in  fact  vanishes  at  the  shadow  boundary.  The  numerical 
variability  (noise)  resulting  from  taking  the  difference  of  two  relatively  large  quantities  of 
vanishing  difference  was  solved  by  examining  the  geometrical  origins  of  the  magnitude 
difference  term  suggested  by  the  equation 

r=p-p'  (17) 

Extracting  square  magnitudes,  one  obtains 

r2  =  p2  +  p'2  -  2/?/?'cos(^-  ft)  (18) 

In  the  vicinity  of  the  shadow  boundary  <f>  -  <j)'  =  n  +  e  and  the  above  equation  may  be 
rewritten  as 


r2  =  p2  +  p'2  -  2/?//cos(;r  +  s)  (19) 

Expanding  the  cosine  function  for  small  displacements  e  about  tc  ,  one  obtains 
-  1  +  s2  /  2  and  the  above  equation  becomes 

r2  =  p2  +  p'2  -  2pp\-  1  +  £2  / 2)  (20) 

r2  =  p2  +  p'2  -  2 pp'  -  pp's2 
r 2  =  [p+  p2)-  pp'e2 


Expanding  the  radical  for  [pp's  2)«(P  +  p1)2 ,  one  obtains 


PP  ? 

r  =  p\  p'  -  — - 


2  {p+  P') 


(21) 


and  the  quantity  of  interest  becomes 


PPL  p 

2  (p+  P') 


p+p'-r  = 


(22) 


and  Eq.  [16]  may  be  written 


Ed  = 


y~ikr j 


where 


dr  =  p\  p'  - 


ppL  p 

2  {p+  ft) 


(23) 


(24) 


Eq.  (23)  is  of  course  exact  and  the  approximation  in  Eq.  (24)  is  only  required  in  the 
vicinity  of  the  shadow  boundary  and  was  found  to  effectively  eliminate  numerical 
variability  (noise)  in  plots  of  the  interfering  fields. 


Figure  36 


13.1.2.2.  Comparison  with  other  models 

A  graphical  comparison  of  path  loss  contour  plots  in  the  vertical  propagation  plane 
obtained  from  four  established  propagation  models  and  of  the  proposed  model  in  its 
current  state  of  development  appears  in  Figure  37  and  Figure  38.  The  most  striking 
difference  among  the  coverage  presentations  in  the  lit  region  is  the  relative  presence  of 
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the  fine  interference  structure  as  seen  in  APM  and  VTRPE.  Though  striking,  such 
differences  are  usually  merely  indicative  of  whether  the  model  in  question  suppresses 
superposition  of  direct  and  ground  specular  reflections.  The  profile  essentially  involves 
three  perfectly  absorbing  knife  edges  in  the  first,  third  and  fourth  positions  which  create 
three  successive  shadows  asymptotically  and  a  fourth  knife  edge  in  the  second  position 
which  is  only  a  partial  mask  for  the  for  the  third  knife  edge  and  casts  only  a  finite 
shadow.  A  vertical  excursion  at  great  distances  for  such  a  configuration  should  encounter 
three  shadow  boundaries  and  should  exhibit  a  monotonic  increase  in  loss  upon  traversing 
the  first  shadow  boundary  from  above.  This  behavior  is  indeed  exhibited  in  APM, 
VTRPE,  and  our  proposed  model,  whereas  ITM  exhibits  precipitous  discontinuities  at  the 
first  and  second  shadow  boundaries  and  fails  to  detect  a  third  shadow  boundary  and  the 
TIREM  exhibits  a  misplaced  region  of  excessive  loss  in  the  vicinity  of  the  second  shadow 
boundary. 

A  horizontal  excursion  proceeding  from  the  transmitter  for  such  a  terrain  profile  is 
expected  to  exhibit  a  monotonic  increase  in  loss  upon  penetrating  the  first  shadow 
boundary  and  to  then  experience  two  additional  shadow  boundaries  if  at  an  elevation 
greater  than  the  third  knife  edge  and  one  additional  shadow  if  at  a  lesser  elevation.  These 
features  are  again  exhibited  in  the  APM,  VTRPE,  and  the  proposed  model  whereas 
TIREM  again  has  difficulties  in  the  second  shadow  for  low  and  high  elevation  horizontal 
excursions  and  ITM  model  exhibits  non-physically  severe  losses  across  both 
interpretations  of  its  first  shadow  which  curiously  intersect  at  the  third  knife  edge. 

Thus  unlike  TIREM  and  ITM,  our  model  agrees  well  with  the  more  rigorous  results 
found  in  APM  and  VTRPE  with  the  advantage  that  it’s  run  time  are  comparable  with 
TIREM  and  ITM. 
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